---
title: Policy Leakage or Policy Benefit? Spatial Spillovers from Conservation Policies
  in Common Property Resources
author: Mani Rouhi Rad, Dale T. Manning, Jordan F. Suter, Chris Goemans
affiliation: Clemson University, Colorado State University
abstract: "Voluntary, incentive-based policies are often offered to resource users to address the over-exploitation of common property resources. Spillovers from these policies have important implications for policy effectiveness. Considering different externalities that can affect producers' irrigation decisions, we investigate whether permanently retiring groundwater rights impacts groundwater use among active neighboring groundwater users. We find that enrollment in the retirement program causes neighboring wells to initially pump less groundwater on average. However, water use reductions are only temporary, due to relative increases in stock levels over time. The results imply that policies that retire water rights may conserve more water than anticipated in the short term, but, over time, higher resource stocks could lead to policy leakage. We also present evidence that the decrease in groundwater use is driven by the shared resource stock. Our results provide empirical evidence for the importance of spatial spillovers related to common property resource management."
date: '`r format(Sys.Date(), "%B %d, %Y")`'
documentclass: article
output: 
  bookdown::pdf_document2:
    keep_tex: yes
    number_sections: yes
bibliography: conservation_leakage_ref.bib 
link-citations: true
linkcolor: blue
fontsize: 12pt
header-includes:
    - \usepackage{setspace}\doublespacing
    - \usepackage{float}
    - \usepackage{lscape}
    - \newcommand{\blandscape}{\begin{landscape}}
    - \newcommand{\elandscape}{\end{landscape}}
    - \usepackage{chngcntr}
    - \counterwithin{figure}{section}
    - \counterwithin{table}{section}

---


```{r global_options, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache=TRUE, fig.pos = 'H')
```


```{r, include=FALSE, warning=FALSE}
library(foreign)
library(data.table)
library(plyr)
library(dplyr)
library(reshape2)
library(AER)
library(ggplot2)
library(ggsn)
library(ggspatial)
library(rgeos)
library(rgdal)
library(raster)
library(scales)
library(sp)
library(sf)
library(maptools)
library(RColorBrewer)
library(broom)
require(stargazer)
library(knitr)
library(plm)
library(lfe)
library(multcomp)
library(prism)
```

```{r, include=FALSE, warning=FALSE}




setwd("/Users/manirouhirad/Dropbox/OWCAP Econ/Code/CREP_paper/Submission code and files")

#~~~~~~~~~~~~~~~~~~~~~~~~~~
# Insert data
#~~~~~~~~~~~~~~~~~~~~~~~~~~

# Read WIMAS data
filenames = list.files(path = "./input_data/WIMAS_csv" , pattern="*.csv", full.names=TRUE)
ldf <- lapply(filenames, read.csv)
ldf <- mapply(cbind, ldf, "year"=filenames, SIMPLIFY=F)
WIMAS_dt = rbindlist(ldf)
WIMAS_dt[, WUR_CODE:= as.character(WUR_CODE)]
WIMAS_dt[, WUR_CODE:= gsub(" ", "", WUR_CODE)]
setkey(WIMAS_dt, WUR_CODE)
WIMAS_dt[, year := NULL]

# Arkansas RB CREP
crep = readOGR(dsn =  "./input_data/Arkansas_RB_CREP", layer = "CREP_bnd_060306_tr")
crep              <- spTransform(crep, CRS("+proj=utm +zone=14 ellps=WGS84")) # transorm to utm for distance units
crep_buffer <- gUnaryUnion(crep)
crep_buffer = gBuffer(crep_buffer, byid = TRUE, width = 8046.72) # 5 miles


# Counties
counties = readOGR("./input_data/UScounties", layer = "UScounties")
counties = counties[counties@data$STATE_NAME %like% "Kansas", ]
counties <- spTransform(counties, "+proj=utm +zone=14 ellps=WGS84")

# High Plains
high_plains <- readOGR(dsn="./input_data/HP_boudary", layer="hp_bound2010")
high_plains <-  spTransform(high_plains, "+proj=utm +zone=14 ellps=WGS84")
high_plains <- raster::crop(high_plains, counties)

# Hydraulic conductivity
HP_K =  readOGR(dsn = "./input_data/HP_hydraulic_conductivity", layer = "HP_K_range")
HP_K <- spTransform(HP_K, "+proj=utm +zone=14 ellps=WGS84")

# specific yield
HP_SY = raster("./input_data/hp_sp_yield/SY_HP.tif")
HP_SY = setMinMax(HP_SY)
HP_SY = projectRaster(HP_SY, crs = "+proj=utm +zone=14 ellps=WGS84")

#========================================================================================
# # # # # # # #                      Data Management                      # # # # # # # #
#========================================================================================

baz= WIMAS_dt[complete.cases(WUA_YEAR) & WUA_YEAR>1995]
baz[, WR_NUM_QUAL := paste(WR_NUM, WR_QUAL, sep = "-")]
baz[, PRIORITY := as.Date(PRIORITY, "%d-%b-%Y")]
baz=unique(baz, by=c("WR_ID", "WRF_ACTIVE", "RIGHT_TYPE", "VCNTY_CODE", "WR_NUM", "WR_QUAL", "UMW_CODE", "WRF_STATUS",
                     "SOURCE", "PRIORITY", "PDIV_ID", "FPV_ACTIVE", "TWP", "TWP_DIR", "RNG", "RNG_DIR",
                     "SECT", "DWR_ID", "FEET_NORTH", "FEET_WEST", "QUAL_FOUR", "QUAL_THREE", "QUAL_TWO", "QUAL_ONE",
                     "FO_NUM", "BASIN_NUM", "BASIN_NAME", "GWMD_NUM", "CNTY_ABREV", "CNTY_NAME", "SUA_NAME", "STREAM_NUM",
                     "NUM_WELLS", "LOT_NUMBER", "LOT_QUAL1", "LOT_QUAL2", "FPDIV_COMM", "WUAFO_NUM", "WUAPERS_ID", "WUACOR_NUM",
                     "WUA_YEAR", "WUAUMWCODE", "HOURS_PUMP", "PUMP_RATE", "METER_QTY", "METER_UNIT", "WUR_CODE", "RPT_WRIGHT",
                     "ACRES_IRR", "DATE_MEAS", "DPTH_WATER", "DPTH_WELL", "REEL_NUM", "BLIP_NUM", "RPT_DATE", "CHEM_IND",
                     "CROP_CODE", "B_METER_RD", "E_METER_RD", "SYSTEM", "AF_USED", "FPDIV_KEY", "WRF_KEY", "LONGITUDE",
                     "LATITUDE", "QUANT_ID", "AUTH_QUANT", "ADD_QUANT", "QUANT_UNIT", "QSTOR_IND", "RATE_ID", "AUTH_RATE",
                     "ADD_RATE", "RATE_UNIT", "RSTOR_IND", "TACRES_IRR", "NACRES_IRR", "LAST_NAME", "FIRST_NAME", "CHG",
                     "FILE_ID"))

baz[, location := paste(LONGITUDE, LATITUDE, TWP, TWP_DIR, RNG, RNG_DIR, SECT, DWR_ID, FEET_NORTH, FEET_WEST, sep = ", ")]

setkey(baz, WR_NUM, PDIV_ID, WUA_YEAR, CNTY_NAME)
baz_all = copy(baz)

# CREP wells: This required a little data cleaning by hand. 
# The basic idea here is that some CREP wells used to have a different PDIV_ID. This part of the code is making that connection.

# new PDIV_ID assigned
# assigned CREP but zero acre-feet in history of the well.
baz[WR_NUM == 1198, PDIV_ID := 67138] # new PDIV_ID assigned
baz[WR_NUM == 1198, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2018 ONLY; SEND WUR THRU 2018"]

baz = baz[WR_NUM != 8884] # assigned CREP but zero acre-feet in history of the well.

baz[WR_NUM == 9496, PDIV_ID := 75566] # new PDIV_ID assigned
baz[WR_NUM == 9496, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2018 ONLY; SEND WUR THRU 2018"]

baz[WR_NUM == 10437, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2020 ONLY - SEND WUR THROUGH 2021"]

baz = baz[WR_NUM!=10521 | (PDIV_ID == 68477 | PDIV_ID == 75544)]

baz = baz[!(WR_NUM == 12631 & WR_ID == 12887 & PDIV_ID == 18601)]
baz = baz[!(WR_NUM == 12651 & WR_ID == 12909 & PDIV_ID == 18601)]

baz[WR_NUM == 9496 & WR_ID == 13054, PDIV_ID := 47973] # new PDIV_ID assigned

baz[WR_NUM == 13540, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2020 ONLY - SEND WUR THROUGH 2021"]

baz[WR_NUM == 16320 & WR_ID == 16634, PDIV_ID := 72731] # new PDIV_ID assigned

baz[WR_NUM == 16320, FPDIV_COMM := "CREP IRR THROUGH 12/31/2010 ONLY"]

baz = baz[PDIV_ID != 14177]
baz = baz[PDIV_ID != 14177]
baz = baz[WR_NUM  != 19973]

baz[WR_NUM == 20820 & WR_ID == 21201, PDIV_ID := 81961] # new PDIV_ID assigned
baz[WR_NUM == 20820 & WR_ID == 21201, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2015 ONLY"]

baz[WR_NUM == 20939 & WR_ID == 21320, PDIV_ID := 83190] # new PDIV_ID assigned
baz[WR_NUM == 20939 & WR_ID == 21320, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2020 ONLY - SEND WUR THROUGH 2021"]

baz[WR_NUM == 20940 & WR_ID == 21321, PDIV_ID := 78654] # new PDIV_ID assigned
baz[WR_NUM == 20940 & WR_ID == 21321, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2020 ONLY - SEND WUR THROUGH 2021"]

baz[WR_NUM == 20941 & WR_ID == 21322, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2020 ONLY; SEND WUR THROUGH 2020"]

baz[WR_NUM == 20942 & WR_ID == 21323, PDIV_ID := 70154] # new PDIV_ID assigned
baz[WR_NUM == 20942 & WR_ID == 21323, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2020 ONLY; SEND WUR THROUGH 2020"]

baz[WR_NUM == 21335 & WR_ID == 21720,  PDIV_ID := 68732] # new PDIV_ID assigned
baz[WR_NUM == 21335 & WR_ID == 21720 , FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010"]

baz = baz[PDIV_ID != 24618]

baz[WR_NUM == 23012 & WR_ID == 23410, PDIV_ID := 62017] # new PDIV_ID assigned

baz[WR_NUM == 23012 & WR_ID == 23410 & PDIV_ID == 62017, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2013 ONLY"]

baz[WR_NUM == 23344 & WR_ID == 23751, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2020 ONLY; SEND WUR THROUGH 2020"]

baz[WR_NUM == 25256 & WR_ID == 25685 & PDIV_ID == 27551, PDIV_ID := 67301] # new PDIV_ID assigned
baz[WR_NUM == 25256 & WR_ID == 25685 & PDIV_ID == 67301, FPDIV_COMM := "CREP IRR THROUGH 12/31/2010 ONLY"]

baz[WR_NUM  == 25632 & WR_ID == 26068 & PDIV_ID == 32065, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2018 ONLY - SEND WUR THRU 2018"]
baz[WR_NUM  == 25632 & WR_ID == 26068 & PDIV_ID == 53032, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2018 ONLY - SEND WUR THRU 2018"]

baz = baz[!(WR_NUM == 27521 & WR_ID == 27979 & PDIV_ID == 23401)]

baz[WR_NUM == 27521 & WR_ID == 27979 & (PDIV_ID == 48152 | PDIV_ID == 77518), PDIV_ID := 48152] # new PDIV_ID assigned
baz[WR_NUM == 27521 & WR_ID == 27979 &  PDIV_ID == 48152, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2015 ONLY"]

baz[WR_NUM == 28283 & WR_ID == 28746 & PDIV_ID == 30947, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2019 ONLY - SENT WUR THRU 2019"]

baz[WR_NUM == 30691 & WR_ID == 31166 & PDIV_ID == 45803, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2020 ONLY; SEND WUR THROUGH 2020"]
baz[WR_NUM == 30692 & WR_ID == 31167 & PDIV_ID == 40201, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2020 ONLY; SEND WUR THROUGH 2020"]

baz[WR_NUM == 33592 & WR_ID == 34112 & PDIV_ID == 50357, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2020 ONLY - SEND WUR THROUGH 2021"]

baz[WR_NUM == 33599 & WR_ID == 34119,  PDIV_ID := 72952]
baz[WR_NUM == 33599 & WR_ID == 34119 & PDIV_ID == 72952, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2013 ONLY"]

baz[WR_NUM == 33612 & WR_ID == 34132 & PDIV_ID == 28735, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2020 ONLY - SEND WUR THROUGH 2021"]

baz = baz[WR_NUM != 20137043]

baz[WR_NUM == 33620 & WR_ID == 34140,  PDIV_ID := 67339] # new PDIV_ID assigned
baz[WR_NUM == 33620 & WR_ID == 34140,  FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010"]

baz[WR_NUM == 33635 & WR_ID == 34155,  PDIV_ID := 67714] # new PDIV_ID assigned
baz[WR_NUM == 33635 & WR_ID == 34155,  FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2018 ONLY; SEND WUR THRU 2018"]

baz = baz[!(WR_NUM == 36722 & WR_ID == 32768 & PDIV_ID == 5171)]
baz = baz[!(WR_NUM == 36722 & WR_ID == 32768 & PDIV_ID == 11055)]
baz = baz[!(WR_NUM == 36722 & WR_ID == 32768 & PDIV_ID == 12006)]
baz = baz[!(WR_NUM == 36722 & WR_ID == 32768 & PDIV_ID == 19459)]
baz = baz[!(WR_NUM == 36722 & WR_ID == 32768 & PDIV_ID == 25942)]
baz = baz[!(WR_NUM == 36722 & WR_ID == 32768 & PDIV_ID == 29643)]
baz = baz[!(WR_NUM == 36722 & WR_ID == 32768 & PDIV_ID == 31722)]
baz = baz[!(WR_NUM == 36722 & WR_ID == 32768 & PDIV_ID == 41069)]

baz[WR_NUM == 37498 & WR_ID == 38055,  FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2018 ONLY - SEND WUR THRU 2018"]
baz[WR_NUM == 41212 & WR_ID == 41854,  FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2019 ONLY - SEND WUR THRU 2019"]

baz = baz[!(WR_NUM == 41212 & WR_ID == 41854 & PDIV_ID == 30947)]

baz = baz[!(WR_NUM == 9 & WR_ID == 50279 & PDIV_ID == 18601)]

baz[WR_NUM == 19726 & WR_ID == 20091 & PDIV_ID == 26025,  WR_ID := 70721]
baz[WR_NUM == 19726 & WR_ID == 20091 & PDIV_ID == 70721, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010 ONLY"]
baz[WR_NUM == 20526 & WR_ID == 20907 & PDIV_ID == 54607,  WR_ID := 70722]
baz[WR_NUM == 20526 & WR_ID == 70722 & PDIV_ID == 54607, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010"]
baz[WR_NUM == 20526 & WR_ID == 20907 & PDIV_ID == 25788,  WR_ID := 70723]
baz[WR_NUM == 20526 & WR_ID == 70723 & PDIV_ID == 25788, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2012 ONLY"]
baz[WR_NUM == 20526 & WR_ID == 20907 & PDIV_ID == 44086,  WR_ID := 70724]
baz[WR_NUM == 20526 & WR_ID == 70724 & PDIV_ID == 44086, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010"]
baz[WR_NUM == 20526 & WR_ID == 20907 & PDIV_ID == 54367,  WR_ID := 70725]
baz[WR_NUM == 20526 & WR_ID == 70725 & PDIV_ID == 54367, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010"]
baz[WR_NUM == 23447 & WR_ID == 23860 & PDIV_ID == 17180,  WR_ID := 70732]
baz[WR_NUM == 23447 & WR_ID == 70732 & PDIV_ID == 17180, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010 ONLY"]
baz[WR_NUM == 23662 & WR_ID == 24078 & PDIV_ID == 14572,  WR_ID := 70734]
baz[WR_NUM == 23662 & WR_ID == 70734 & PDIV_ID == 14572, FPDIV_COMM := "CREP IRRIGAITON THROUGH 12/31/2010 ONLY"]
baz[WR_NUM == 25383 & WR_ID == 25813 & PDIV_ID == 38636,  WR_ID := 70737]
baz[WR_NUM == 25383 & WR_ID == 70737 & PDIV_ID == 38636, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010 ONLY"]
baz[WR_NUM == 25383 & WR_ID == 25813 & PDIV_ID == 6005,  WR_ID := 70740]
baz[WR_NUM == 25383 & WR_ID == 70740 & PDIV_ID == 6005, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010 ONLY"]
baz[WR_NUM == 25383 & WR_ID == 25813 & PDIV_ID == 72363,  WR_ID := 70743]
baz[WR_NUM == 25383 & WR_ID == 25813 & PDIV_ID == 16208,  WR_ID := 70743]
baz[WR_NUM == 25383 & WR_ID == 70743 & PDIV_ID == 16208,  PDIV_ID := 72363]
baz[WR_NUM == 25383 & WR_ID == 70743 & PDIV_ID == 72363, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010 ONLY"]
baz[WR_NUM == 25386 & WR_ID == 25816 & PDIV_ID == 41060,  WR_ID := 70749]
baz[WR_NUM == 25386 & WR_ID == 70749 & PDIV_ID == 41060, FPDIV_COMM := "CREP IRRIGAITON THROUGH 12/31/2010"]
baz[WR_NUM == 25386 & WR_ID == 25816 & PDIV_ID == 4379,  WR_ID := 70750]
baz[WR_NUM == 25386 & WR_ID == 70750 & PDIV_ID == 4379, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010 ONLY"]
baz[WR_NUM == 29170 & WR_ID == 29639 & PDIV_ID == 21813,  WR_ID := 70751]
baz[WR_NUM == 29170 & WR_ID == 70751 & PDIV_ID == 21813, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010"]
baz[WR_NUM == 25385 & WR_ID == 25815 & PDIV_ID == 33105,  WR_ID := 70755]
baz[WR_NUM == 25385 & WR_ID == 70755 & PDIV_ID == 33105, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010 ONLY"]
baz[WR_NUM == 25385 & WR_ID == 25815 & PDIV_ID == 24623,  PDIV_ID := 72381]
baz[WR_NUM == 25385 & WR_ID == 25815 & PDIV_ID == 72381,  WR_ID := 70756]
baz[WR_NUM == 25385 & WR_ID == 70756 & PDIV_ID == 72381, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010 ONLY"]
baz[WR_NUM == 25385 & WR_ID == 25815 & PDIV_ID == 11976,  WR_ID := 70758]
baz[WR_NUM == 25385 & WR_ID == 70758 & PDIV_ID == 11976, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010 ONLY"]
baz[WR_NUM == 25385 & WR_ID == 25815 & PDIV_ID == 46472,  WR_ID := 70759]
baz[WR_NUM == 25385 & WR_ID == 70759 & PDIV_ID == 46472, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010 ONLY"]
baz[WR_NUM == 29170 & WR_ID == 29639 & PDIV_ID == 40026,  WR_ID := 70763]
baz[WR_NUM == 29170 & WR_ID == 70763 & PDIV_ID == 40026, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010 ONLY"]
baz[WR_NUM == 19088 & WR_ID == 19446 & PDIV_ID == 30623,  WR_ID := 70950]
baz[WR_NUM == 19088 & WR_ID == 70950 & PDIV_ID == 30623, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010"]

baz = baz[!(WR_NUM == 19088 & WR_ID == 70951 & PDIV_ID == 30623)]

baz[WR_NUM == 27342 & WR_ID == 27800 & PDIV_ID == 19673,  WR_ID := 70957]
baz[WR_NUM == 27342 & WR_ID == 70957 & PDIV_ID == 19673, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010 ONLY"]
baz[WR_NUM == 23724 & WR_ID == 24141 & PDIV_ID == 52306,  WR_ID := 71064]
baz[WR_NUM == 23724 & WR_ID == 71064 & PDIV_ID == 52306, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2010 ONLY"]
baz[WR_NUM == 28307 & WR_ID == 28774 & PDIV_ID == 26996,  WR_ID := 71394]
baz[WR_NUM == 28307 & WR_ID == 71394 & PDIV_ID == 26996, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2011 ONLY"]
baz[WR_NUM == 24790 & WR_ID == 27948 & PDIV_ID == 65492,  WR_ID := 71435]
baz[WR_NUM == 24790 & WR_ID == 71435 & PDIV_ID == 65492, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2011 ONLY"]
baz[WR_NUM == 24790 & WR_ID == 27948 & PDIV_ID == 3606,  WR_ID := 71436]
baz[WR_NUM == 24790 & WR_ID == 71436 & PDIV_ID == 3606, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2011 ONLY"]
baz[WR_NUM == 18273 & WR_ID == 18613 & PDIV_ID == 22870,  WR_ID := 71291]
baz[WR_NUM == 18273 & WR_ID == 71291 & PDIV_ID == 22870, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2011 ONLY"]
baz[WR_NUM == 20804 & WR_ID == 21185 & PDIV_ID == 42709,  WR_ID := 72291]
baz[WR_NUM == 20804 & WR_ID == 72291 & PDIV_ID == 42709, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2012 ONLY"]
baz[WR_NUM == 21833 & WR_ID == 22222 & PDIV_ID == 41067,  WR_ID := 72370]
baz[WR_NUM == 21833 & WR_ID == 72370 & PDIV_ID == 41067, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2011 ONLY"]
baz[WR_NUM == 25967 & WR_ID == 26404 & PDIV_ID == 53633,  WR_ID := 72913]
baz[WR_NUM == 25967 & WR_ID == 72913 & PDIV_ID == 53633, FPDIV_COMM := "CREP IRRIGTION THROUGH 12/31/2018 ONLY (SEND WUR THROUGH 2018)"]
baz[WR_NUM == 25967 & WR_ID == 26404 & PDIV_ID == 18956,  WR_ID := 72915]
baz[WR_NUM == 25967 & WR_ID == 72915 & PDIV_ID == 18956, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2020 ONLY - SEND WUR THRU 2020"]
baz[WR_NUM == 25967 & WR_ID == 26404 & PDIV_ID == 17109,  WR_ID := 72916]
baz[WR_NUM == 25967 & WR_ID == 72916 & PDIV_ID == 17109, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2020 ONLY - SEND WUR THRU 2020"]
baz[WR_NUM == 25968 & WR_ID == 26405 & PDIV_ID == 27054,  WR_ID := 72920]
baz[WR_NUM == 25968 & WR_ID == 72920 & PDIV_ID == 27054, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2020 ONLY - SEND WUR THRU 2020"]
baz[WR_NUM == 25968 & WR_ID == 26405 & PDIV_ID == 72514,  WR_ID := 72921]
baz[WR_NUM == 25968 & WR_ID == 72921 & PDIV_ID == 72514, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2020 ONLY - SEND WUR THRU 2020"]
baz[WR_NUM == 24174 & WR_ID == 24953 & PDIV_ID == 10076,  WR_ID := 75965]
baz[WR_NUM == 24174 & WR_ID == 74965 & PDIV_ID == 10076, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2013 ONLY"]
baz[WR_NUM == 18293 & WR_ID == 18633 & PDIV_ID == 51290,  WR_ID := 18633]
baz[WR_NUM == 18293 & WR_ID == 75314 & PDIV_ID == 51290, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2014 ONLY"]
baz[WR_NUM == 21384 & WR_ID == 21769 & PDIV_ID == 2463,  WR_ID := 75445]
baz[WR_NUM == 21384 & WR_ID == 75445 & PDIV_ID == 2463, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2017 ONLY"]
baz[WR_NUM == 21384 & WR_ID == 21769 & PDIV_ID == 11228,  WR_ID := 75480]
baz[WR_NUM == 21384 & WR_ID == 75480 & PDIV_ID == 11228, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2017 ONLY"]
baz[WR_NUM == 21384 & WR_ID == 21769 & PDIV_ID == 16686,  WR_ID := 75485]
baz[WR_NUM == 21384 & WR_ID == 75485 & PDIV_ID == 16686, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2017 ONLY"]
baz[WR_NUM == 21384 & WR_ID == 75487 & PDIV_ID == 26488,  PDIV_ID := 80307]
baz[WR_NUM == 21384 & WR_ID == 21769 & PDIV_ID == 26488,  PDIV_ID := 80307]
baz[WR_NUM == 21384 & WR_ID == 21769 & PDIV_ID == 80307,  WR_ID   := 75487]
baz[WR_NUM == 21384 & WR_ID == 75487 & PDIV_ID == 80307,  FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2017 ONLY"]
baz[WR_NUM == 21384 & WR_ID == 21769 & PDIV_ID == 35508,  WR_ID := 75488]
baz[WR_NUM == 21384 & WR_ID == 75488 & PDIV_ID == 35508, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2017 ONLY"]
baz[WR_NUM == 14275 & WR_ID == 14554 & PDIV_ID == 46462,  WR_ID := 75502]
baz[WR_NUM == 14275 & WR_ID == 75502 & PDIV_ID == 46462, FPDIV_COMM := "CREP IRRIGATION THROUGH 12-31-2013 ONLY"]
baz[WR_NUM == 19156 & WR_ID == 19514 & PDIV_ID == 33625,  WR_ID := 76990]
baz[WR_NUM == 19156 & WR_ID == 76990 & PDIV_ID == 33625, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2013 ONLY"]
baz[WR_NUM == 18072 & WR_ID == 18411 & PDIV_ID == 75672,  WR_ID := 78271]
baz = baz[!(WR_NUM == 18072 & WR_ID == 18411 & PDIV_ID == 66277) | WUA_YEAR < 2008]
baz[WR_NUM == 18072 & WR_ID == 18411 & PDIV_ID == 66277,  PDIV_ID := 75672]
baz[WR_NUM == 18072 & WR_ID == 18411 & PDIV_ID == 75672,  WR_ID := 78271]
baz[WR_NUM == 18072 & WR_ID == 78271 & PDIV_ID == 75672,  FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2014 ONLY"]
baz[WR_NUM == 23662 & WR_ID == 24078 & PDIV_ID == 45206,  WR_ID := 78562]
baz[WR_NUM == 23662 & WR_ID == 70733 & PDIV_ID == 45206,  WR_ID := 78562]
baz[WR_NUM == 23662 & WR_ID == 78562 & PDIV_ID == 45206, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2014 ONLY"]
baz[WR_NUM == 23662 & WR_ID == 24078 & PDIV_ID == 38575,  WR_ID := 82573]
baz[WR_NUM == 23662 & WR_ID == 70733 & PDIV_ID == 38575,  WR_ID := 82573]
baz[WR_NUM == 23662 & WR_ID == 82573 & PDIV_ID == 38575, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2017 ONLY; SEND WUR THROUGH 2018"]
baz[WR_NUM == 26067 & WR_ID == 26504 & PDIV_ID == 44203,  WR_ID := 83436]
baz[WR_NUM == 26067 & WR_ID == 83436 & PDIV_ID == 44203, FPDIV_COMM := "CREP IRRIGATION THROUGH 12/31/2018 ONLY (SEND WUR THROUGH 2018)"]
baz[WR_NUM == 22575 & WR_ID == 22965 & PDIV_ID == 54868,  WR_ID := 83747]
baz[WR_NUM == 22575 & WR_ID == 83747 & PDIV_ID == 54868, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2020 ONLY; SEND WUR THROUGH 2020"]
baz[WR_NUM == 21056 & WR_ID == 83783 & PDIV_ID == 49449,  WR_ID := 83783]
baz[WR_NUM == 21056 & WR_ID == 21436 & PDIV_ID == 49449, FPDIV_COMM := "CREP IRRIGATION THRU 12/31/2018 ONLY; SEND WUR THRU 2018"]


# create a dataset of CREP wells
baz_CREP_WIMAS = baz[FPDIV_COMM %like% "CREP"]
baz_CREP_WIMAS[, signed_year := FPDIV_COMM]
baz_CREP_WIMAS[FPDIV_COMM %like% "THROUGH" | FPDIV_COMM %like% "THRU"]
baz_CREP_WIMAS[, foo_gh := regexpr("THROUGH", signed_year)]
baz_CREP_WIMAS[, foo_u  := regexpr("THRU",    signed_year)]
baz_CREP_WIMAS[foo_gh >0 & foo_u >0]
baz_CREP_WIMAS[, foo := ifelse(foo_u<0, foo_gh, foo_u)]
baz_CREP_WIMAS[, signed_year := as.character(signed_year)]
baz_CREP_WIMAS[, signed_year := substr(signed_year, foo, nchar(signed_year))]
baz_CREP_WIMAS[, signed_year := gsub("THROUGH ", "", signed_year)]
baz_CREP_WIMAS[, signed_year := gsub("THRU ", "", signed_year)]
baz_CREP_WIMAS[, signed_year := substr(signed_year, 1, 10)]
baz_CREP_WIMAS[, signed_year := as.Date(signed_year, "%m/%d/%Y")]

setkey(baz_CREP_WIMAS, WR_ID, PDIV_ID, WUA_YEAR)
baz_CREP_WIMAS[, group:= .GRP, by=c("WR_ID", "PDIV_ID")]
baz_CREP_WIMAS[, c("FPDIV_COMM", "foo_gh", "foo_u", "foo") := NULL]

baz_CREP_WIMAS[, foo := max(WUA_YEAR), by="group"]
baz_CREP_WIMAS[, long_foo := ifelse(WUA_YEAR == foo, LONGITUDE, 0)]
baz_CREP_WIMAS[, lat_foo  := ifelse(WUA_YEAR == foo, LATITUDE,  0)]

baz_CREP_WIMAS[, long_foo := min(long_foo), by="group"]
baz_CREP_WIMAS[, lat_foo  := max(lat_foo) , by="group"]

baz_CREP_WIMAS[, LONGITUDE := long_foo]
baz_CREP_WIMAS[, LATITUDE  := lat_foo]
baz_CREP_WIMAS[, c("foo", "long_foo", "lat_foo") := NULL]


# take care of the GEO CTR: this process basically considers if water use is zero for the wells and non-zero for the GEO CTR, then allocates GEO CTR for each well and sets GEO_CENTERED == 1 so that
# we can unique these wells to 1 well.
baz[, PDIV_ID_GC    := ifelse(FPDIV_COMM %like% "GEO CTR", PDIV_ID, 0)]
baz[, PDIV_ID_GC    := max(PDIV_ID_GC), by=c("WR_ID")]
baz[, water_use_GC  := ifelse(PDIV_ID == PDIV_ID_GC, AF_USED, 0)]
baz[, water_use_NGC := ifelse(PDIV_ID != PDIV_ID_GC, AF_USED, 0)]
baz[, water_use_GC  := max(water_use_GC), by=c("WR_ID", "WUA_YEAR")]

baz[, water_use     := ifelse(water_use_GC > 0 & water_use_NGC == 0, water_use_GC, AF_USED)]
baz[, GEO_CENTERED  := ifelse(water_use_GC > 0 & water_use_NGC == 0, 1, 0)]
baz[, GEO_CENTERED  := max(GEO_CENTERED), by=c("WR_ID", "WUA_YEAR")]
baz = baz[GEO_CENTERED == 0 | FPDIV_COMM %like% "GEO CTR"] # this removes the non-geo centers.

# remove wells that are always zero ...
baz[, foo2 := sum(water_use), by=c("WR_ID", "PDIV_ID")]
baz = baz[foo2>0]

baz_CREP_WIMAS[, last_year := ifelse(AF_USED>0, 1, 0)] # this can be water_use > 10 for example
baz_CREP_WIMAS[, last_year_2 := cumsum(last_year), by=c("WR_ID", "PDIV_ID")]
baz_CREP_WIMAS[  last_year==0, last_year_2 := 0]
baz_CREP_WIMAS[, last_year_3 := max(last_year_2), by=c("WR_ID", "PDIV_ID")]
baz_CREP_WIMAS[, signed_year_2 := ifelse(last_year_2==last_year_3, WUA_YEAR, 0)]
baz_CREP_WIMAS[, signed_year_2 := max(signed_year_2), by=c("WR_ID", "PDIV_ID")]
baz_CREP_WIMAS[, c("last_year", "last_year_2", "last_year_3", "group") := NULL]

# signed_year_2 is the last year that a well extracted groundwater.

## Remove the CREP wells from baz, the dataset of active wells
foo = baz_CREP_WIMAS[,.(WR_NUM, WR_ID, PDIV_ID, id=1)]
foo = unique(foo)

setkey(foo, WR_NUM, WR_ID, PDIV_ID)
setkey(baz, WR_NUM, WR_ID, PDIV_ID)

baz = foo[baz]
baz = baz[is.na(id)]
baz[, id := NULL]

# wells within CREP well buffers
baz_positive = baz[water_use>0]
baz_all_unique_sp = unique(baz_positive, by=c("WR_ID", "PDIV_ID"))
coord             = baz_all_unique_sp[,.(LONGITUDE, LATITUDE)]
baz_all_unique_sp = baz_all_unique_sp[,.(WR_ID, PDIV_ID)]
baz_all_unique_sp = SpatialPointsDataFrame(coord, data = baz_all_unique_sp, # I got the idea here: https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf
                                           proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
baz_all_unique_sp <- spTransform(baz_all_unique_sp, CRS("+proj=utm +zone=14 ellps=WGS84")) # transorm to utm for distance units
crep_buffer       <- spTransform(crep_buffer, CRS("+proj=utm +zone=14 ellps=WGS84")) # transorm to utm for distance units
crep              <- spTransform(crep, CRS("+proj=utm +zone=14 ellps=WGS84")) # transorm to utm for distance units


## active wells within the UARB buffer area. This is the study area
wells_within_ARB    <- over(baz_all_unique_sp,    crep_buffer)
wells_within_ARB    <- spCbind(baz_all_unique_sp, wells_within_ARB)
wells_within_ARB_data= data.table(coordinates(wells_within_ARB), wells_within_ARB@data)
wells_within_ARB_data= wells_within_ARB_data[complete.cases(wells_within_ARB)]
wells_within_ARB_data= wells_within_ARB_data[,.(WR_ID, PDIV_ID)]


setkey(wells_within_ARB_data, WR_ID, PDIV_ID)
setkey(baz_positive,          WR_ID, PDIV_ID)
wells_within_ARB_data = baz_positive[wells_within_ARB_data]


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Histogram of CREP wells over time

# To generate this figure, you need to run the code above until line with: wells_within_ARB_data = baz_positive[wells_within_ARB_data]
# except that line with crep_buffer = gBuffer(crep_buffer, byid = TRUE, width = 8046.72) # 5 miles
# needs to change to 0 width of buffer:
# crep_buffer = gBuffer(crep_buffer, byid = TRUE, width = 0)
# and also run the following:

# foo = unique(wells_within_ARB_data[AF_USED > 1 & ACRES_IRR > 1], by=c("LONGITUDE", "LATITUDE"))
# foo = nrow(foo)
# baz_CREP_sp = baz_CREP_WIMAS[,.(WR_ID, PDIV_ID, WUA_YEAR,
#                                 WUAPERS_ID, LAST_NAME, FIRST_NAME, AF_USED, ACRES_IRR, signed_year, signed_year_2, LONGITUDE, LATITUDE)]
# baz_CREP_sp = unique(baz_CREP_sp, by=c("WR_ID", "PDIV_ID"))
# coord = baz_CREP_sp[,.(LONGITUDE, LATITUDE)]
# baz_CREP_sp = SpatialPointsDataFrame(coord, data = baz_CREP_sp,
#                                      proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
# baz_CREP_sp <- spTransform(baz_CREP_sp, CRS("+proj=utm +zone=14 ellps=WGS84")) # transorm to utm for distance units
# foo3_app=data.table(group= as.numeric(as.character(rownames(baz_CREP_sp@data))), signed_year= baz_CREP_sp@data$signed_year, signed_year_2= baz_CREP_sp@data$signed_year_2)
# 
# dt_1 = foo3_app[, .N, by="signed_year_2"]
# setkey(dt_1, signed_year_2)
# dt_1[, N := cumsum(N)]
# 
# dt_1[, ratio := round(N/foo*100, 2)]
# dt_1 = dt_1[signed_year_2 > 2007]
# dt_1 = melt(dt_1, "signed_year_2")
# 
# dt_1[variable=="ratio", value := value * 20]
# dt_1[variable== "N", variable := "Cumulative"]
# dt_1[variable== "ratio", variable := "Percent"]

dt_1 = fread("./regression_data/data_for_CREP_histogram.csv")


fig_hist = ggplot(dt_1, aes(x=as.factor(signed_year_2), y=value, group=variable, fill=variable))+
  geom_col(position = "dodge", width = .7)+
  theme_bw()+
  scale_y_continuous(sec.axis = sec_axis(~.*.05, name = "Percent of wells in CREP"))+
  scale_fill_manual("", values=c("grey20", "grey70"))+
  xlab("Year")+
  ylab("Number of wells in CREP")+
  theme(axis.text    = element_text(size=13),
        axis.title   = element_text(size=16),
        axis.text.x  = element_text(angle = 90, hjust = 1),
        legend.position   = "bottom",
        legend.text       = element_text(size=16),
        legend.title      = element_text(size=16),
        legend.key.size   = unit(.7,"cm")
  )  
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


coord = wells_within_ARB_data[,.(LONGITUDE, LATITUDE)]
wells_within_ARB_sp = SpatialPointsDataFrame(coord, data = wells_within_ARB_data[,.(WR_ID, PDIV_ID, WUA_YEAR)], # I got the idea here: https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf
                                             proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
wells_within_ARB_sp <- spTransform(wells_within_ARB_sp, CRS("+proj=utm +zone=14 ellps=WGS84")) # transorm to utm for distance units



#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# CREP wells and CREP area
baz_CREP_sp_fig = unique(baz_CREP_WIMAS[signed_year_2>2007], by=c("WR_ID", "PDIV_ID", "LONGITUDE", "LATITUDE"))
coord       = baz_CREP_sp_fig[,.(LONGITUDE, LATITUDE)]
baz_CREP_sp_fig = baz_CREP_sp_fig[,.(WR_ID, PDIV_ID)]
baz_CREP_sp_fig = SpatialPointsDataFrame(coord, data = baz_CREP_sp_fig, # I got the idea here: https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf
                                         proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
baz_CREP_sp_fig <- spTransform(baz_CREP_sp_fig, CRS("+proj=utm +zone=14 ellps=WGS84")) # transorm to utm for distance units

UARB_river = readOGR("./input_data/Arkansas_River", layer = "Arkansas_River_KS")
UARB_river = spTransform(UARB_river, "+proj=utm +zone=14 ellps=WGS84")
UARB_river <- raster::crop(UARB_river, crep_buffer)


UARB_river_sf       = st_as_sf(UARB_river)
crep_buffer_sf      = st_as_sf(crep_buffer)
wells_within_ARB_sf = st_as_sf(wells_within_ARB_sp)
baz_CREP_sf         = st_as_sf(baz_CREP_sp_fig)


fig_CREP_wells = ggplot() +
  geom_sf(data = wells_within_ARB_sf, col = "grey80") +
  geom_sf(data = baz_CREP_sf, col="black", size=3) +
  geom_sf(data = crep_buffer_sf, fill=NA, size=1, col="black", linetype ="dashed") +
  geom_sf(data = UARB_river_sf, col="grey30", size=1) +
  north(crep_buffer_sf, location = "topleft", scale = .2, symbol = 12) +
  annotation_scale(location = "bl", width_hint = 0.4) +
  theme(
    axis.ticks.y      = element_blank(),
    axis.text.y       = element_blank(),
    axis.ticks.x      = element_blank(),
    axis.text.x       = element_blank(),
    axis.title        = element_blank(),
    panel.background  = element_rect(fill='white'))


high_plains   = st_as_sf(high_plains)
counties      = st_as_sf(counties)

fig_CREP_area  = ggplot() +
  geom_sf(data = counties,      fill=NA, col = "black") +
  geom_sf(data = high_plains,   fill="grey80") +
  geom_sf(data = crep_buffer_sf,fill="grey30", size=1, col="black") +
  theme(
    axis.ticks.y      = element_blank(),
    axis.text.y       = element_blank(),
    axis.ticks.x      = element_blank(),
    axis.text.x       = element_blank(),
    axis.title        = element_blank(),
    panel.background  = element_rect(fill='white'))

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Merge with Saturated Thickness
coord=wells_within_ARB_data[,.(LONGITUDE, LATITUDE)]
wells_within_ARB_data_sp = SpatialPointsDataFrame(coord, data = wells_within_ARB_data,
                                                  proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
wells_within_ARB_data_sp = spTransform(wells_within_ARB_data_sp, "+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-101 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")


filenames=list.files(path = "input_data/HPA_sat_thickness_inputs" , pattern="*.tif", full.names=TRUE)
filenames=filenames[seq(1, length(filenames), by=3)]

foo_4=data.table()
for (i in 1:length(filenames)) {
  assign(paste0("sat_thickness", 1999+i), raster(filenames[i]))
  foo      = get(paste0("sat_thickness", 1999+i))
  foo[foo < 0] <- 0
  foo      = foo * 3.28084
  foo_2    = wells_within_ARB_data_sp[wells_within_ARB_data_sp@data$WUA_YEAR == 2000+i, ]
  foo_3    = data.frame(foo_2@data,
                        raster::extract(foo, foo_2))
  foo_3=data.table(foo_3)
  foo_4=rbind(foo_4, foo_3)
  # print(i)
}

foo_4_saturated_thickness=foo_4[raster..extract.foo..foo_2.>=0]


# Merge with PRISM data

# GDD                     

foo_4_saturated_thickness[, id_temp := 1:.N]

# This part is commented out because it will take time to run, but in order to generate GDD, the following code was used here:

# library(prism)
# foo_4_saturated_thickness_2 = foo_4_saturated_thickness[WUA_YEAR>2000]
# foo_4_saturated_thickness_2[, id_temp := 1:.N]
# 
# options(prism.path = "OWCAP Econ/Data/PRISM/daily/temp/mean/")
# get_prism_dailys(type="tmean", minDate = "1996-04-01", maxDate = "2017-10-01", keepZip=F)
# prism_path=ls_prism_data(absPath = TRUE)
# prism_path=data.table(prism_path)
# prism_path=prism_path[,.(abs_path)]
# # prism_path=prism_path[c(1:95, 97:105)]
# 
# prism_path[, foo := regexpr("\\_[^\\_]*$", abs_path)]
# prism_path[, baz := substr(abs_path, 1, foo-1)]
# prism_path[, foo := regexpr("\\_[^\\_]*$", baz)]
# prism_path[, baz := substr(abs_path, foo+1, nchar(baz))]
# prism_path[, month := as.numeric(substr(baz, 5, 6))]
# prism_path[, day   := as.numeric(substr(baz, 7, 8))]
# prism_path[, year  := as.numeric(substr(baz, 1, 4))]
# 
# prism_path=prism_path[month > 3 & month < 10 & year>2000]
# prism_path=prism_path[,.(abs_path)]
# 
# a=data.table()
# for (i in 1:(nrow(prism_path))) {
#   path_r=prism_path[i]
#   raster_prism=raster(path_r$abs_path)
#   
#   g <- regexpr("\\_[^\\_]*$", path_r$abs_path)
#   g <- substr(path_r$abs_path, 1, g[1]-1)
#   h <- regexpr("\\_[^\\_]*$", g)
#   h <- substr(g, h[1]+1, nchar(g))
#   g <- as.numeric(substr(h, 5, 6))
#   k <- as.numeric(substr(h, 7, 8))
#   h <- as.numeric(substr(h, 1, 4))
#   
#   panel_saturated_thickness_foo = foo_4_saturated_thickness_2[WUA_YEAR==h, .(id_temp, LONGITUDE, LATITUDE)]
#   coord=panel_saturated_thickness_foo[,.(LONGITUDE, LATITUDE)]
#   panel_saturated_thickness_foo = SpatialPointsDataFrame(coord, panel_saturated_thickness_foo[,.(id_temp)], 
#                                                          proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
#   panel_saturated_thickness_foo <- spTransform(panel_saturated_thickness_foo, crs(raster_prism))
#   
#   t_mean = data.table(raster::extract(raster_prism, panel_saturated_thickness_foo))
#   panel_saturated_thickness_foo = data.table(panel_saturated_thickness_foo@data, t_mean)
#   panel_saturated_thickness_foo[, day := k]
#   panel_saturated_thickness_foo[, month := g]
#   panel_saturated_thickness_foo[, year := h]
#   a=rbind(a, panel_saturated_thickness_foo)
#   print(i)
# }
# 
# 
# setkey(a, id_temp, year, month, day)
# a[, V1 := ifelse(V1 < 10, 10, V1)]
# a[, V1 := ifelse(V1 > 30, 30, V1)]
# a = a[month < 10]
# a[, V1 := V1 - 10]
# a[, GDD := sum(V1), by=c("id_temp", "year")]
# a = unique(a, by=c("id_temp", "year"))
# GDD_data = a[,.(id_temp, GDD)]
# 
# saveRDS(GDD_data, "OWCAP Econ/Data/PRISM/GDD_data_RR.rds")


GDD_data = readRDS("input_data/PRISM/GDD_data_RR.rds")

setkey(foo_4_saturated_thickness, id_temp)
setkey(GDD_data, id_temp)
foo_4_saturated_thickness= foo_4_saturated_thickness[GDD_data]
a = foo_4_saturated_thickness[,.(LONGITUDE, LATITUDE, id_temp, WUA_YEAR)]

# Precipitation

options(prism.path = "input_data/PRISM/precip/")
get_prism_monthlys(type="ppt", year = 2000:2017, mon = 4:9, keepZip=F)
prism_path=ls_prism_data(absPath = TRUE)
prism_path=data.table(prism_path)
prism_path=prism_path[,.(abs_path)]

b=data.table()
for (i in 1:(nrow(prism_path))) {
  path_r=prism_path[i]
  raster_prism=raster(path_r$abs_path)
  
  g <- regexpr("\\_[^\\_]*$", path_r$abs_path)
  g <- substr(path_r$abs_path, 1, g[1]-1)
  h <- regexpr("\\_[^\\_]*$", g)
  h <- substr(g, h[1]+1, nchar(g))
  g <- as.numeric(substr(h, 5, 6))
  h <- as.numeric(substr(h, 1, 4))
  
  panel_saturated_thickness_foo = a[WUA_YEAR==h]
  coord=panel_saturated_thickness_foo[,.(LONGITUDE, LATITUDE)]
  panel_saturated_thickness_foo = SpatialPointsDataFrame(coord, panel_saturated_thickness_foo, 
                                                         proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
  panel_saturated_thickness_foo <- spTransform(panel_saturated_thickness_foo, crs(raster_prism))
  
  precip = data.table(raster::extract(raster_prism, panel_saturated_thickness_foo))
  panel_saturated_thickness_foo = data.table(panel_saturated_thickness_foo@data, precip)
  b=rbind(b, panel_saturated_thickness_foo)
  # print(i)
}

setnames(b, old = "V1", new = "precip")

b[, precip := mean(precip), by=c("id_temp", "WUA_YEAR")]
b= unique(b, by=c("id_temp", "WUA_YEAR"))


coord=b[,.(LONGITUDE, LATITUDE)]
saturated_thickness_temp_precip = b[,.(id_temp, precip)]
setkey(GDD_data, id_temp)
setkey(saturated_thickness_temp_precip, id_temp)
saturated_thickness_temp_precip = saturated_thickness_temp_precip[GDD_data]

foo_4_saturated_thickness[, saturated_thickness := raster..extract.foo..foo_2.]
foo_4_saturated_thickness[, c("raster..extract.foo..foo_2.", "GDD") := NULL]

setkey(saturated_thickness_temp_precip, id_temp)
setkey(foo_4_saturated_thickness, id_temp)

wells_within_buffer_radius_panel_2 = foo_4_saturated_thickness[saturated_thickness_temp_precip]
wells_within_buffer_radius_panel_2 = wells_within_buffer_radius_panel_2[complete.cases(saturated_thickness)]

# Add hydraulic conductivity, K, and specific yield, SY.
coord = wells_within_buffer_radius_panel_2[,.(LONGITUDE, LATITUDE)]
wells_within_buffer_radius_panel_sp = SpatialPointsDataFrame(coord, data = wells_within_buffer_radius_panel_2, 
                                                             proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
wells_within_buffer_radius_panel_sp <- spTransform(wells_within_buffer_radius_panel_sp, CRS("+proj=utm +zone=14 ellps=WGS84")) # transorm to utm for distance units

# Hydraulic conductivity

b_sp_K=over(wells_within_buffer_radius_panel_sp, HP_K)
b_sp_K= data.table(b_sp_K)

wells_within_buffer_radius_panel_2 = data.table(wells_within_buffer_radius_panel_2, b_sp_K[,.(range_K=RANGE, lower_K=MAJOR1, upper_K=MINOR1)])

# specific yield

wells_within_buffer_radius_panel_2    = data.table(wells_within_buffer_radius_panel_2,
                                                   sp_yield = raster::extract(HP_SY, wells_within_buffer_radius_panel_sp))

# find the distance between CREP wells and active wells
# basically, this part finds info such as number of wells, water use, acres, etc for all the wells that are not in CREP but are in the UARB region that are in the
# XX mile radius of our treated wells. The idea is to control for their decisions in our regressions

baz_CREP_WIMAS = baz_CREP_WIMAS[signed_year_2 > 2007]

baz_CREP_sp = baz_CREP_WIMAS[,.(WR_ID, PDIV_ID, WUA_YEAR,
                                WUAPERS_ID, LAST_NAME, FIRST_NAME, AF_USED, ACRES_IRR, signed_year, signed_year_2, LONGITUDE, LATITUDE, AF_USED)]
baz_CREP_sp[, AF_USED := max(AF_USED), by=c("LONGITUDE", "LATITUDE")]
baz_CREP_sp = unique(baz_CREP_sp, by=c("WR_ID", "PDIV_ID"))

coord = baz_CREP_sp[,.(LONGITUDE, LATITUDE)]
baz_CREP_sp = SpatialPointsDataFrame(coord, data = baz_CREP_sp,
                                     proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
baz_CREP_sp <- spTransform(baz_CREP_sp, CRS("+proj=utm +zone=14 ellps=WGS84")) # transorm to utm for distance units

wells_within_buffer_radius_panel_2[, id := 1:.N]
coord   = wells_within_buffer_radius_panel_2[,.(LONGITUDE, LATITUDE)]
foo1_sp = SpatialPointsDataFrame(coord, data = wells_within_buffer_radius_panel_2[,.(id)],
                                 proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
foo1_sp <- spTransform(foo1_sp, CRS("+proj=utm +zone=14 ellps=WGS84")) # transorm to utm for distance units
# foo1_sp = gBuffer(foo1_sp, byid = TRUE, width = 8046.72*1/5) # 1 mile

foo3=data.table()
i=1
for (i in 1:nrow(data.table(baz_CREP_sp@data))) {
  foo2_sp = baz_CREP_sp[i,]
  foo2= data.table(distance = spDistsN1(foo1_sp, foo2_sp, longlat = FALSE))
  foo2[, id := 1:.N]
  foo3 = rbind(foo3, foo2)
  # print(i)
}

foo3[, group:= ifelse(id==1, 1, 0)]
foo3[, group:= cumsum(group)]
foo3[, distance := distance/(8046.72*1/5)]

foo3_app=data.table(group= as.numeric(as.character(rownames(baz_CREP_sp@data))), signed_year= baz_CREP_sp@data$signed_year, signed_year_2= baz_CREP_sp@data$signed_year_2, AF_USED = baz_CREP_sp@data$AF_USED, CREP_Person_ID = baz_CREP_sp@data$WUAPERS_ID)

setkey(foo3, group)
setkey(foo3_app, group)

CREP_distance = foo3[foo3_app]
setkey(CREP_distance, id)
rm(foo3_app, foo3)

foo = CREP_distance[, min(distance), by="id"]

foo[, dist_bin := ifelse(V1 == 0, 0,
                         ifelse(V1 <= 1, 1,
                                ifelse(V1 <= 2, 2,
                                       ifelse(V1 <= 3, 3,
                                              ifelse(V1 <= 4, 4, 5)))))]

foo = foo[,.(id, dist_bin)]
setkey(foo, id)
setkey(wells_within_buffer_radius_panel_2, id)
wells_within_buffer_radius_panel_2 = wells_within_buffer_radius_panel_2[foo]
wells_within_buffer_radius_panel_2 = wells_within_buffer_radius_panel_2[dist_bin != 0]
wells_within_buffer_radius_panel_2[, well_id := paste(WR_ID, PDIV_ID, sep = "-")]
wells_within_buffer_radius_panel_2[, county_year := paste(CNTY_ABREV, WUA_YEAR, sep = "-")]
# wells_within_buffer_radius_panel_2[water_use>0]

wells_within_buffer_radius_panel_2_b = copy(wells_within_buffer_radius_panel_2)

foo = baz_CREP_WIMAS[ ,.(WUAPERS_ID, id=1)]
foo = unique(foo)
setkey(foo, WUAPERS_ID)
setkey(wells_within_buffer_radius_panel_2_b, WUAPERS_ID)

wells_within_buffer_radius_panel_2_b = foo[wells_within_buffer_radius_panel_2_b]
wells_within_buffer_radius_panel_2_b = wells_within_buffer_radius_panel_2_b[is.na(id)]

wells_within_buffer_radius_panel_2_b[, id := i.id]
wells_within_buffer_radius_panel_2_b[, i.id := NULL]

rm(list = c('sat_thickness1995','sat_thickness1996','sat_thickness1997','sat_thickness1998','sat_thickness1999','sat_thickness2000'
            ,'sat_thickness2001','sat_thickness2002','sat_thickness2003','sat_thickness2004','sat_thickness2005','sat_thickness2006'
            ,'sat_thickness2007','sat_thickness2008','sat_thickness2009','sat_thickness2010','sat_thickness2011','sat_thickness2012'
            ,'sat_thickness2013','sat_thickness2014','sat_thickness2015','sat_thickness2016'))

rm(list = c('b_HP','b_hp_sp','b_HP_U','b_sp','b_sp_K','b2'
            ,'baseline_temp_precip','baz','baz_all','baz_all_unique_sp','baz_positive'
            ,'c','coord','coord_b','coord_b_hp','coord_CO','coord_KS'
            ,'coord_merge','coord_NE','coord_OK','counties','counties_for_coord'))

rm(list = c('foo_4_saturated_thickness','foo1_sp','high_plains','foo2_sp', 'saturated_thickness_temp_precip', 't_mean'
            ,'ldf', 'a', 'b',  'foo', 'foo_2', 'foo_3', 'foo_4', 
            'panel_saturated_thickness_foo', 'path_r', 'precip', 'wells_within_ARB', 'wells_within_ARB_data', 'wells_within_ARB_data_sp', 'wells_within_ARB_sp'))


wells_within_buffer_radius_panel_2[,   range_year   := paste(TWP, TWP_DIR, RNG, RNG_DIR, WUA_YEAR, sep = "-")]
wells_within_buffer_radius_panel_2_b[, range_year   := paste(TWP, TWP_DIR, RNG, RNG_DIR, WUA_YEAR, sep = "-")]

wells_within_buffer_radius_panel_2_c = copy(wells_within_buffer_radius_panel_2_b)
wells_within_buffer_radius_panel_2_c[, section_year := paste(TWP, TWP_DIR, RNG, RNG_DIR, SECT, WUA_YEAR, sep = "-")]
wells_within_buffer_radius_panel_2_c[, section_only := paste(TWP, TWP_DIR, RNG, RNG_DIR, SECT, sep = "-")]
wells_within_buffer_radius_panel_2_c = wells_within_buffer_radius_panel_2_c[WUA_YEAR>2000]

# saveRDS(wells_within_buffer_radius_panel_2_b, "regression_data/wells_within_buffer_radius_panel_2_b.rds")

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Motivation figure

foo = wells_within_buffer_radius_panel_2_b[water_use > 0,.(dist_bin, WUA_YEAR, water_use)]
foo[, dist_bin_2 := ifelse(dist_bin <= 2 , "< 2 mi" , 
                           ifelse(dist_bin <= 4 , "2-4 mi" , "> 4 mi"))]
foo = foo[, mean(water_use), by=c("dist_bin_2", "WUA_YEAR")]
foo[, dist_bin_2 := factor(dist_bin_2, ordered = T, levels = c("< 2 mi", "2-4 mi", "> 4 mi"))]

fig_motivation = ggplot(foo, aes(x=WUA_YEAR, y=V1, group= factor(dist_bin_2), col= factor(dist_bin_2) ))+
  geom_point(size=2)+
  geom_line(size=2)+
  geom_vline(xintercept = 2008) + 
  theme_classic()+
  xlab("Year")+
  ylab("Water use (AF)")+
  scale_color_brewer()+
  theme(axis.text    = element_text(size=12),
        axis.title   = element_text(size=12),
        legend.position   = "bottom",
        legend.text       = element_text(size=12),
        legend.title      = element_text(size=12),
        legend.key.size   = unit(1,"cm")  )  +
  guides(col=guide_legend(title="Distance (mile)"))


# Summary statistics
a= copy(wells_within_buffer_radius_panel_2_b)
a = unique(a, by=c("LONGITUDE", "LATITUDE"))
a = a[,.N, by="dist_bin"]

with_2mi    = a[dist_bin == 1 | dist_bin ==2, sum(N)] # wells inside 2 mile buffer
outside_2mi = a[dist_bin == 3 | dist_bin ==4 | dist_bin ==5, sum(N)] # wells outside 2 mile buffer

# 773+591+20 # this is for the hydraulic conductivity for comparison

a= wells_within_buffer_radius_panel_2_b[water_use>0 & WUA_YEAR < 2008,.(dist_bin, water_use, saturated_thickness, precip, GDD, lower_K, upper_K, ACRES_IRR)]
a[, ACRES_IRR := as.numeric(ACRES_IRR)]
a[, lower_K := as.numeric(as.character(lower_K))]
a[, upper_K := as.numeric(as.character(upper_K))]
a[, dist_bin := ifelse(dist_bin < 3, 1 , 2)]
cols = colnames(a)[2:ncol(a)]

a_count = a[,.N, by="dist_bin"]
a_count[, N := sqrt(N)]

b = a[, lapply(.SD, function(x) sd(x, na.rm=T)), by="dist_bin", .SDcols = cols]
a[, (cols) := lapply(.SD, function(x) mean(x, na.rm=T)), by="dist_bin", .SDcols = cols]
a = unique(a, by="dist_bin")
setkey(b,       dist_bin)
setkey(a_count, dist_bin)
b = b[a_count]

cols = colnames(a)[2:ncol(a)]
b[, (cols) := lapply(.SD, "/", sqrt(N)), by="dist_bin", .SDcols = cols]

count_foo = wells_within_buffer_radius_panel_2_b[water_use>0 & WUA_YEAR == 2007, .N, by="dist_bin"]
count_foo = count_foo[dist_bin == 1 | dist_bin == 2]
b[, N := NULL]
setkey(b,         dist_bin)
setkey(count_foo, dist_bin)
b = b[count_foo]


a[, (cols) := lapply(.SD, function(x) round(x, 2)), .SDcols = cols]
b[, (cols) := lapply(.SD, function(x) round(x, 2)), .SDcols = cols]
cols = data.table(colnames(b))
cols[, V2 := paste(V1, "SD", sep = "_")]
setnames(b, old = cols$V1, new = cols$V2)

setkey(a, dist_bin)
setkey(b, dist_bin_SD)
a = a[b]

a= data.table(melt(a, "dist_bin"))
sum_stats = reshape(a, idvar = "variable", timevar = "dist_bin", direction = "wide")


sum_stats_2 = readRDS("regression_data/CREP_wells_sum_stats.rds") # this is done for CREP wells similar to active wells.
count_CREP_wells = nrow(unique(sum_stats_2, by=c("PDIV_ID", "WR_ID")))

sum_stats_2 = sum_stats_2[AF_USED > 0 & WUA_YEAR < 2008,.(water_use = AF_USED, saturated_thickness, precip, GDD, lower_K, upper_K, ACRES_IRR)]
sum_stats_2[, ACRES_IRR := as.numeric(ACRES_IRR)]
sum_stats_2[, lower_K := as.numeric(as.character(lower_K))]
sum_stats_2[, upper_K := as.numeric(as.character(upper_K))]
cols = colnames(sum_stats_2)

a_count = sqrt(sum_stats_2[,.N])
b = sum_stats_2[, lapply(.SD, function(x) sd(x, na.rm=T)), .SDcols = cols]
cols = colnames(b)
b = data.table(id=1, b[, lapply(.SD, "/", a_count), .SDcols = cols])
b[, N:= count_CREP_wells]


sum_stats_2[, (cols) := lapply(.SD, function(x) mean(x, na.rm=T)), .SDcols = cols]
sum_stats_2 = unique(sum_stats_2)
sum_stats_2 = data.table(id=1, sum_stats_2)

sum_stats_2[, (cols) := lapply(.SD, function(x) round(x, 2)), .SDcols = cols]
b[, (cols) := lapply(.SD, function(x) round(x, 2)), .SDcols = cols]
cols = data.table(colnames(b))
cols[, V2 := paste(V1, "SD", sep = "_")]
setnames(b, old = cols$V1, new = cols$V2)

setkey(sum_stats_2, id)
setkey(b, id_SD)
sum_stats_2 = sum_stats_2[b]

sum_stats_2 = data.table(melt(sum_stats_2, "id"))
sum_stats_2[, id := NULL]

sum_stats[variable == "N_SD", value.1 := with_2mi]
sum_stats[variable == "N_SD", value.2 := outside_2mi]

setkey(sum_stats,   variable)
setkey(sum_stats_2, variable)

sum_stats = sum_stats_2[sum_stats]
setnames(sum_stats, old = c("value", "value.1", "value.2"), new = c("CREP wells", "< 2 miles", "> 2 miles"))

sum_stats[variable == "water_use",           variable := "Groundwater extraction (AF)"]
sum_stats[variable == "saturated_thickness", variable := "Saturated thickness (ft)"]
sum_stats[variable == "precip",              variable := "Precipitation (mm)"]
sum_stats[variable == "GDD",                 variable := "Growing degree days"]
sum_stats[variable == "lower_K",             variable := "Hydraulic conductivity (lower bound)"]
sum_stats[variable == "upper_K",             variable := "Hydraulic conductivity (upper bound)"]
sum_stats[variable == "ACRES_IRR",           variable := "Acres irrigated"]
sum_stats[, id := c(1, 3, 5, 7, 9, 11, 13, 2, 4, 6, 8, 10, 12, 14, 15)]
setkey(sum_stats, id)
sum_stats[, `CREP wells` := ifelse(id== 2 | id== 4 | id== 6 | id== 8 | id== 10 | id== 12 | id== 14, paste0("(", `CREP wells`, ")"), `CREP wells`)]
sum_stats[, `< 2 miles`  := ifelse(id== 2 | id== 4 | id== 6 | id== 8 | id== 10 | id== 12 | id== 14, paste0("(", `< 2 miles`,  ")"), `< 2 miles` )]
sum_stats[, `> 2 miles`  := ifelse(id== 2 | id== 4 | id== 6 | id== 8 | id== 10 | id== 12 | id== 14, paste0("(", `> 2 miles`,  ")"), `> 2 miles` )]
sum_stats[variable == "N_SD",           variable := "Number of wells"]
sum_stats[variable %like% "_SD", variable := ""]
sum_stats[, id := NULL]

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Person_ID_CREP = CREP_distance[distance <= 2 ,.(id, CREP_Person_ID)]
Person_ID_CREP = unique(Person_ID_CREP)
Person_ID_CREP = Person_ID_CREP[, lapply(.SD, paste0, collapse=","), by = id]


wells_within_buffer_radius_panel_2_all = copy(wells_within_buffer_radius_panel_2)


wells_within_buffer_radius_panel_2_b = wells_within_buffer_radius_panel_2_b[,.(WR_ID, WR_NUM, PDIV_ID, LONGITUDE, LATITUDE, WUAPERS_ID, WUA_YEAR)]
wells_within_buffer_radius_panel_2_b[, qq:=1]
# all

setkey(wells_within_buffer_radius_panel_2,   WR_ID, WR_NUM, PDIV_ID, LONGITUDE, LATITUDE, WUAPERS_ID, WUA_YEAR)
setkey(wells_within_buffer_radius_panel_2_b, WR_ID, WR_NUM, PDIV_ID, LONGITUDE, LATITUDE, WUAPERS_ID, WUA_YEAR)

wells_within_buffer_radius_panel_2 = merge(wells_within_buffer_radius_panel_2_b, wells_within_buffer_radius_panel_2, 
                                           by=c("WR_ID", "WR_NUM", "PDIV_ID", "LONGITUDE", "LATITUDE", "WUAPERS_ID", "WUA_YEAR"), all = T)

wells_within_buffer_radius_panel_2 = wells_within_buffer_radius_panel_2[is.na(qq)]

foo = data.table()
qux = data.table()
wells_within_buffer_radius_panel_3 = data.table()

#---
wells_within_buffer_radius_panel_3 = copy(wells_within_buffer_radius_panel_2[WUA_YEAR>2000])
setkey(wells_within_buffer_radius_panel_3, id)

qux = wells_within_buffer_radius_panel_3[,.(id, WUA_YEAR)]
foo = CREP_distance[distance <= 2 & distance > 0]
setkey(foo, id)
setkey(qux, id)
foo= foo[qux]
foo= foo[complete.cases(group)]
foo[, distance := round(distance, 3)]
foo[, count:= ifelse(signed_year_2+1 <= WUA_YEAR, 1, 0)]
foo[, min_distance  := min(distance), by=c("id", "WUA_YEAR")]
foo[, count := sum(count), by=c("id", "WUA_YEAR")]
foo[, AF_USED_2 := ifelse(signed_year_2+1 <= WUA_YEAR, AF_USED, 0)]
foo[, AF_USED_2 := sum(AF_USED_2), by=c("id", "WUA_YEAR")]

foo[AF_USED_2>0]
foo[count>0]

foo[, AF_USED_2 := sum(AF_USED), by=c("id", "WUA_YEAR")]
foo = foo[distance == min_distance]
foo = unique(foo, by=c("id", "WUA_YEAR"))
foo = foo[,.(id, group, distance, count, AF_USED_2)]

setkey(foo, id)
setkey(wells_within_buffer_radius_panel_2, id)
wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_2[foo]


wells_within_buffer_radius_panel_3[, group_year := paste(group, WUA_YEAR, sep = "-")]

wells_within_buffer_radius_panel_3[, long_lat := paste(LONGITUDE, LATITUDE, sep = "-")]
wells_within_buffer_radius_panel_3 = unique(wells_within_buffer_radius_panel_3, by=c("long_lat", "WUA_YEAR"))

wells_within_buffer_radius_panel_3[, year_2 := ifelse(count>0, 1, 0)]



setkey(wells_within_buffer_radius_panel_3, id)
setkey(Person_ID_CREP, id)
wells_within_buffer_radius_panel_3 = Person_ID_CREP[wells_within_buffer_radius_panel_3]
wells_within_buffer_radius_panel_3[, ww := "0"]

i=1
for (i in 1:nrow(wells_within_buffer_radius_panel_3)) {
  WUAPERS_ID_1     = wells_within_buffer_radius_panel_3[i, WUAPERS_ID]
  CREP_Person_ID_1 = wells_within_buffer_radius_panel_3[i, CREP_Person_ID]
  a = grepl(pattern = WUAPERS_ID_1, CREP_Person_ID_1, fixed = T)
  wells_within_buffer_radius_panel_3[i, ww := a]
  # print(i)
}

wells_within_buffer_radius_panel_3[, vv := ifelse(ww==T, 1, 0)]
wells_within_buffer_radius_panel_3[, vv := max(vv), by=c("LONGITUDE", "LATITUDE")]

# saveRDS(wells_within_buffer_radius_panel_3, "regression_data/appendix_reg_table_456.rds")


panel_FE_02_avg_2_org <- felm(water_use ~ year_2 +  GDD + precip | county_year + long_lat | 0 | long_lat,
                              data = wells_within_buffer_radius_panel_3)

panel_FE_02_avg_3_org <- felm(water_use ~ year_2 +  GDD + precip | group_year + long_lat | 0 | long_lat,
                              data = wells_within_buffer_radius_panel_3)

wells_within_buffer_radius_panel_3[count ==0, AF_USED_2 := 0]

panel_FE_02_N_1_2_org <- felm(water_use ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                              data = wells_within_buffer_radius_panel_3)

panel_FE_02_N_1_3_org <- felm(water_use ~ count +  GDD + precip | group_year + long_lat | 0 | long_lat, 
                              data = wells_within_buffer_radius_panel_3)


panel_FE_02_avg_2_1 <- felm(water_use ~ year_2 +  GDD + precip | county_year + long_lat | 0 | long_lat,
                            data = wells_within_buffer_radius_panel_3[vv==1])

panel_FE_02_avg_3_1 <- felm(water_use ~ year_2 +  GDD + precip | group_year + long_lat | 0 | long_lat,
                            data = wells_within_buffer_radius_panel_3[vv==1])

wells_within_buffer_radius_panel_3[count ==0, AF_USED_2 := 0]

panel_FE_02_N_1_2_1 <- felm(water_use ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                            data = wells_within_buffer_radius_panel_3[vv==1])

panel_FE_02_N_1_3_1 <- felm(water_use ~ count +  GDD + precip | group_year + long_lat | 0 | long_lat, 
                            data = wells_within_buffer_radius_panel_3[vv==1])




panel_FE_02_avg_2_0 <- felm(water_use ~ year_2 +  GDD + precip | county_year + long_lat | 0 | long_lat,
                            data = wells_within_buffer_radius_panel_3[vv==0])

panel_FE_02_avg_3_0 <- felm(water_use ~ year_2 +  GDD + precip | group_year + long_lat | 0 | long_lat,
                            data = wells_within_buffer_radius_panel_3[vv==0])

panel_FE_02_N_1_2_0 <- felm(water_use ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                            data = wells_within_buffer_radius_panel_3[vv==0])

panel_FE_02_N_1_3_0 <- felm(water_use ~ count +  GDD + precip | group_year + long_lat | 0 | long_lat, 
                            data = wells_within_buffer_radius_panel_3[vv==0])





# #========================================================================================
# # # # # # # # #                       Main Regressions                    # # # # # # # #
# #========================================================================================

# Add section by year and township-range by year
wells_within_buffer_radius_panel_2 = copy(wells_within_buffer_radius_panel_2_c) # sorry this is a bit confusing. But I switched back to the data without _b.
# wells_within_buffer_radius_panel_2[, id := 1:.N]

coord             = wells_within_buffer_radius_panel_2[,.(LONGITUDE, LATITUDE)]
foo = wells_within_buffer_radius_panel_2[,.(id)]
foo = SpatialPointsDataFrame(coord, data = foo, proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
foo <- spTransform(foo, CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0")) # transorm to utm for distance units

KS_PLSS = readOGR(dsn = "input_data/PLSS/cadastral_PLSSTOWN_ks_3664917_01/cadastral", layer = "plsstownship_a_ks") # Kansas township range data
KS_PLSS <- spTransform(KS_PLSS, CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0")) # transorm to utm for distance units

KS_PLSS_over = over(foo, KS_PLSS)
KS_PLSS_over = data.table(KS_PLSS_over)
KS_PLSS_over = KS_PLSS_over[,.(town, range)]

foo = data.table(foo@data, KS_PLSS_over)

setkey(wells_within_buffer_radius_panel_2, id)
setkey(foo, id)
wells_within_buffer_radius_panel_2 = wells_within_buffer_radius_panel_2[foo]
wells_within_buffer_radius_panel_2[, TR_year := paste(town, range, WUA_YEAR, sep = "-")]


coord             = wells_within_buffer_radius_panel_2[,.(LONGITUDE, LATITUDE)]
foo = wells_within_buffer_radius_panel_2[,.(id)]
foo = SpatialPointsDataFrame(coord, data = foo, proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
foo <- spTransform(foo, CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0")) # transorm to utm for distance units

KS_PLSS = readOGR(dsn = "input_data/PLSS/cadastral_PLSSSECT_ks_3665155_01/cadastral", layer = "plsssection_a_ks")
KS_PLSS <- spTransform(KS_PLSS, CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0")) # transorm to utm for distance units


KS_PLSS_over = over(foo, KS_PLSS)
KS_PLSS_over = data.table(KS_PLSS_over)
KS_PLSS_over = KS_PLSS_over[,.(mtrs)]

foo = data.table(foo@data, KS_PLSS_over)

setkey(wells_within_buffer_radius_panel_2, id)
setkey(foo, id)
wells_within_buffer_radius_panel_2 = wells_within_buffer_radius_panel_2[foo]
wells_within_buffer_radius_panel_2[, section_year_2 := paste(mtrs, WUA_YEAR, sep = "-")]


#________________________________
##   Marginal Effects *2* mile 
#________________________________
coord   = wells_within_buffer_radius_panel_2[,.(LONGITUDE, LATITUDE)]
foo1_sp = SpatialPointsDataFrame(coord, data = wells_within_buffer_radius_panel_2[,.(id_1 = id, WUA_YEAR)],
                                 proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
foo1_sp <- spTransform(foo1_sp, CRS("+proj=utm +zone=14 ellps=WGS84")) # transorm to utm for distance units

coord   = wells_within_buffer_radius_panel_2[WUA_YEAR == 2005,.(LONGITUDE, LATITUDE)]
foo3_sp = SpatialPointsDataFrame(coord, data = wells_within_buffer_radius_panel_2[WUA_YEAR == 2005,.(id_3 = id)],
                                 proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
foo3_sp <- spTransform(foo3_sp, CRS("+proj=utm +zone=14 ellps=WGS84")) # transorm to utm for distance units


foo3=data.table()
i=1
for (i in 1:nrow(data.table(foo3_sp@data))) {
  foo = foo3_sp[i,]
  foo2= data.table(distance = spDistsN1(foo1_sp, foo, longlat = FALSE))
  foo2[, distance := distance/(8046.72*1/5)]
  foo2[, id_1 := 1:.N]
  foo2[, id_3 := foo@data$id_3]
  foo2 = foo2[distance <= 2 & distance > 0]
  foo3 = rbind(foo3, foo2)
  # print(i)
}

foo3[, N_2005 := .N, by="id_1"]
foo3 = unique(foo3,  by="id_1")
number_of_wells_2005 = foo3[,.(id_1, N_2005)]



coord   = wells_within_buffer_radius_panel_2[WUA_YEAR == 2004,.(LONGITUDE, LATITUDE)]
foo3_sp = SpatialPointsDataFrame(coord, data = wells_within_buffer_radius_panel_2[WUA_YEAR == 2004,.(id_3 = id)],
                                 proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
foo3_sp <- spTransform(foo3_sp, CRS("+proj=utm +zone=14 ellps=WGS84")) # transorm to utm for distance units


foo3=data.table()
i=1
for (i in 1:nrow(data.table(foo3_sp@data))) {
  foo = foo3_sp[i,]
  foo2= data.table(distance = spDistsN1(foo1_sp, foo, longlat = FALSE))
  foo2[, distance := distance/(8046.72*1/5)]
  foo2[, id_1 := 1:.N]
  foo2[, id_3 := foo@data$id_3]
  foo2 = foo2[distance <= 2 & distance > 0]
  foo3 = rbind(foo3, foo2)
  # print(i)
}

foo3[, N_2004 := .N, by="id_1"]
foo3 = unique(foo3,  by="id_1")
number_of_wells_2004 = foo3[,.(id_1, N_2004)]

coord   = wells_within_buffer_radius_panel_2[WUA_YEAR == 2003,.(LONGITUDE, LATITUDE)]
foo3_sp = SpatialPointsDataFrame(coord, data = wells_within_buffer_radius_panel_2[WUA_YEAR == 2003,.(id_3 = id)],
                                 proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
foo3_sp <- spTransform(foo3_sp, CRS("+proj=utm +zone=14 ellps=WGS84")) # transorm to utm for distance units


foo3=data.table()
i=1
for (i in 1:nrow(data.table(foo3_sp@data))) {
  foo = foo3_sp[i,]
  foo2= data.table(distance = spDistsN1(foo1_sp, foo, longlat = FALSE))
  foo2[, distance := distance/(8046.72*1/5)]
  foo2[, id_1 := 1:.N]
  foo2[, id_3 := foo@data$id_3]
  foo2 = foo2[distance <= 2 & distance > 0]
  foo3 = rbind(foo3, foo2)
  # print(i)
}

foo3[, N_2003 := .N, by="id_1"]
foo3 = unique(foo3,  by="id_1")
number_of_wells_2003 = foo3[,.(id_1, N_2003)]



coord   = wells_within_buffer_radius_panel_2[WUA_YEAR == 2002,.(LONGITUDE, LATITUDE)]
foo3_sp = SpatialPointsDataFrame(coord, data = wells_within_buffer_radius_panel_2[WUA_YEAR == 2002,.(id_3 = id)],
                                 proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
foo3_sp <- spTransform(foo3_sp, CRS("+proj=utm +zone=14 ellps=WGS84")) # transorm to utm for distance units


foo3=data.table()
i=1
for (i in 1:nrow(data.table(foo3_sp@data))) {
  foo = foo3_sp[i,]
  foo2= data.table(distance = spDistsN1(foo1_sp, foo, longlat = FALSE))
  foo2[, distance := distance/(8046.72*1/5)]
  foo2[, id_1 := 1:.N]
  foo2[, id_3 := foo@data$id_3]
  foo2 = foo2[distance <= 2 & distance > 0]
  foo3 = rbind(foo3, foo2)
  # print(i)
}

foo3[, N_2002 := .N, by="id_1"]
foo3 = unique(foo3,  by="id_1")
number_of_wells_2002 = foo3[,.(id_1, N_2002)]


coord   = wells_within_buffer_radius_panel_2[WUA_YEAR == 2001,.(LONGITUDE, LATITUDE)]
foo3_sp = SpatialPointsDataFrame(coord, data = wells_within_buffer_radius_panel_2[WUA_YEAR == 2001,.(id_3 = id)],
                                 proj4string = CRS("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs"))
foo3_sp <- spTransform(foo3_sp, CRS("+proj=utm +zone=14 ellps=WGS84")) # transorm to utm for distance units


foo3=data.table()
i=1
for (i in 1:nrow(data.table(foo3_sp@data))) {
  foo = foo3_sp[i,]
  foo2= data.table(distance = spDistsN1(foo1_sp, foo, longlat = FALSE))
  foo2[, distance := distance/(8046.72*1/5)]
  foo2[, id_1 := 1:.N]
  foo2[, id_3 := foo@data$id_3]
  foo2 = foo2[distance <= 2 & distance > 0]
  foo3 = rbind(foo3, foo2)
  # print(i)
}

foo3[, N_2001 := .N, by="id_1"]
foo3 = unique(foo3,  by="id_1")
number_of_wells_2001 = foo3[,.(id_1, N_2001)]

setkey(number_of_wells_2001, id_1)
setkey(number_of_wells_2002, id_1)
setkey(number_of_wells_2003, id_1)
setkey(number_of_wells_2004, id_1)
setkey(number_of_wells_2005, id_1)

number_of_wells_2001 = number_of_wells_2001[number_of_wells_2002]
number_of_wells_2001 = number_of_wells_2001[number_of_wells_2003]
number_of_wells_2001 = number_of_wells_2001[number_of_wells_2004]
number_of_wells_2001 = number_of_wells_2001[number_of_wells_2005]

number_of_wells_2001[, N_count := N_2001]
number_of_wells_2001[N_2002 > N_count, N_count := N_2002]
number_of_wells_2001[N_2003 > N_count, N_count := N_2003]
number_of_wells_2001[N_2004 > N_count, N_count := N_2004]
number_of_wells_2001[N_2005 > N_count, N_count := N_2005]

number_of_wells_2001 = number_of_wells_2001[,.(id_1, N_count)]


# -------------------
foo = data.table()
qux = data.table()
wells_within_buffer_radius_panel_3 = data.table()

#-----------------
wells_within_buffer_radius_panel_3 = copy(wells_within_buffer_radius_panel_2[WUA_YEAR>2000])
setkey(wells_within_buffer_radius_panel_3, id)

qux = wells_within_buffer_radius_panel_3[,.(id, WUA_YEAR)]
foo = CREP_distance[distance <= 2 & distance > 0]
setkey(foo, id)
setkey(qux, id)
foo= foo[qux]
foo= foo[complete.cases(group)]
foo[, distance := round(distance, 3)]
foo[, count:= ifelse(signed_year_2+1 <= WUA_YEAR, 1, 0)]
foo[, min_distance  := min(distance), by=c("id", "WUA_YEAR")]
foo[, count := sum(count), by=c("id", "WUA_YEAR")]
foo[, AF_USED_2 := ifelse(signed_year_2+1 <= WUA_YEAR, AF_USED, 0)]
foo[, AF_USED_2 := sum(AF_USED_2), by=c("id", "WUA_YEAR")]

foo[, AF_USED_2 := sum(AF_USED), by=c("id", "WUA_YEAR")]
foo = foo[distance == min_distance]
foo = unique(foo, by=c("id", "WUA_YEAR"))

foo = foo[,.(id, group, distance, count, AF_USED_2)]

setkey(foo, id)
setkey(wells_within_buffer_radius_panel_2, id)
wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_2[foo]


wells_within_buffer_radius_panel_3[, group_year := paste(group, WUA_YEAR, sep = "-")]

wells_within_buffer_radius_panel_3[, long_lat := paste(LONGITUDE, LATITUDE, sep = "-")]
wells_within_buffer_radius_panel_3 = unique(wells_within_buffer_radius_panel_3, by=c("long_lat", "WUA_YEAR"))


wells_within_buffer_radius_panel_3[, ACRES_IRR := as.numeric(ACRES_IRR)]
wells_within_buffer_radius_panel_3[, alfalfa := ifelse(CROP_CODE==1, 1, 0)]
wells_within_buffer_radius_panel_3[, corn    := ifelse(CROP_CODE==2, 1, 0)]
wells_within_buffer_radius_panel_3[, sorghum := ifelse(CROP_CODE==3, 1, 0)]
wells_within_buffer_radius_panel_3[, soybean := ifelse(CROP_CODE==4, 1, 0)]
wells_within_buffer_radius_panel_3[, wheat   := ifelse(CROP_CODE==5, 1, 0)]

wells_within_buffer_radius_panel_3[  !(range_K == "25 to 50" | range_K == "50 to 100" | range_K == "100 to 200") , range_K := NA]
wells_within_buffer_radius_panel_3[, conductivity := ifelse(is.na(range_K), NA, 
                                                            ifelse(range_K=="25 to 50" | range_K=="50 to 100", "low", "high"))]
wells_within_buffer_radius_panel_3[, conductivity :=  factor(conductivity)]
wells_within_buffer_radius_panel_3[, conductivity := relevel(conductivity, ref = "low")]


wells_within_buffer_radius_panel_3[, year_2 := ifelse(count>0, 1, 0)]



setkey(wells_within_buffer_radius_panel_3, id)
setkey(number_of_wells_2005, id_1)

wells_within_buffer_radius_panel_3 = number_of_wells_2001[wells_within_buffer_radius_panel_3]
wells_within_buffer_radius_panel_3[is.na(N_count)]
wells_within_buffer_radius_panel_3[is.na(N_count), N_count := 0]
wells_within_buffer_radius_panel_3[, N_1 := max(N_count), by=c("LONGITUDE", "LATITUDE")]
wells_within_buffer_radius_panel_3[, N_2 := max(N_count), by=c("PDIV_ID")]
wells_within_buffer_radius_panel_3[, N_1 := count/N_1]
wells_within_buffer_radius_panel_3[, N_2 := count/N_2]

# saveRDS(wells_within_buffer_radius_panel_3, "regression_data/main_and_appendix_reg_tables_1.rds")


panel_FE_02_avg <- felm(water_use ~ year_2 +  GDD + precip + saturated_thickness | county_year + long_lat | 0 | long_lat,
                        data = wells_within_buffer_radius_panel_3)

panel_FE_02_avg_2 <- felm(water_use ~ year_2 +  GDD + precip | county_year + long_lat | 0 | long_lat,
                          data = wells_within_buffer_radius_panel_3)

panel_FE_02_avg_3 <- felm(water_use ~ year_2 +  GDD + precip | group_year + long_lat | 0 | long_lat,
                          data = wells_within_buffer_radius_panel_3)


panel_FE_02_avg_5 <- felm(water_use ~ year_2 +  GDD + precip                            | county_year + well_id | 0 | well_id, 
                          data = wells_within_buffer_radius_panel_3)

wells_within_buffer_radius_panel_3[count ==0, AF_USED_2 := 0]
panel_FE_02_avg_6 <- felm(water_use ~ AF_USED_2 +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                          data = wells_within_buffer_radius_panel_3)



panel_FE_02_N_1_2 <- felm(water_use ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                          data = wells_within_buffer_radius_panel_3)

panel_FE_02_N_1_3 <- felm(water_use ~ count +  GDD + precip | group_year + long_lat | 0 | long_lat, 
                          data = wells_within_buffer_radius_panel_3)

panel_FE_02_N_1_5 <- felm(water_use ~ count +  GDD + precip                            | county_year + well_id | 0 | well_id, 
                          data = wells_within_buffer_radius_panel_3)

panel_FE_02_N_1_6 <- felm(water_use ~ N_1 +  GDD + precip                              | county_year + PDIV_ID | 0 | PDIV_ID, 
                          data = wells_within_buffer_radius_panel_3[N_1 < 5])



wells_within_buffer_radius_panel_3[, count_2 := ifelse(count==0, 0, 
                                                       ifelse(count==1, 1, 
                                                              ifelse(count==2, 2, 
                                                                     ifelse(count==3, 3, 
                                                                            ifelse(count==4, 4, 
                                                                                   ifelse(count==5, 5, 6))))))]



panel_FE_02_N_1_2_factor <- felm(water_use ~ factor(count_2) +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                                 data = wells_within_buffer_radius_panel_3)




# saturated thickness
panel_FE_02_N_1_ST <- felm(saturated_thickness ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                           data = wells_within_buffer_radius_panel_3)


# # acres irrigated
panel_FE_02_N_1_ACRES_IRR <- felm(ACRES_IRR ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                                  data = wells_within_buffer_radius_panel_3)


panel_FE_02_N_1_ACRES_IRR_2 <- felm(water_use ~ count +  GDD + precip + ACRES_IRR | county_year + long_lat | 0 | long_lat, 
                                    data = wells_within_buffer_radius_panel_3)


# alfalfa
panel_FE_02_N_1_alfalfa <- felm(alfalfa ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                                data = wells_within_buffer_radius_panel_3)

# corn
panel_FE_02_N_1_corn <- felm(corn ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                             data = wells_within_buffer_radius_panel_3)

# sorghum 
panel_FE_02_N_1_sorghum <- felm(sorghum ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                                data = wells_within_buffer_radius_panel_3)

# soybean
panel_FE_02_N_1_soybean <- felm(soybean ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                                data = wells_within_buffer_radius_panel_3)

# wheat
panel_FE_02_N_1_wheat = felm(wheat ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                             data = wells_within_buffer_radius_panel_3)


# hydraulic conductivity lower
wells_within_buffer_radius_panel_3[, .N, by="conductivity"]

panel_FE_02_cond = felm(water_use ~ count +  GDD + precip + conductivity:count | county_year + long_lat | 0 | long_lat, 
                        data = wells_within_buffer_radius_panel_3)

panel_FE_02_cond_3 = felm(water_use ~ year_2 +  GDD + precip + conductivity:year_2 | county_year + long_lat | 0 | long_lat, 
                          data = wells_within_buffer_radius_panel_3)

#### Until 2016 with saturated thickness

wells_within_buffer_radius_panel_4 = readRDS("regression_data/wells_within_buffer_radius_panel_MS_2016.rds")

panel_FE_02_avg_2016 <- felm(water_use ~ year_2 +  GDD + precip + saturated_thickness | county_year + long_lat | 0 | long_lat,
                             data = wells_within_buffer_radius_panel_4)

panel_FE_02_avg_2_2016 <- felm(water_use ~ year_2 +  GDD + precip | county_year + long_lat | 0 | long_lat,
                               data = wells_within_buffer_radius_panel_4)

panel_FE_02_N_1_2016 <- felm(water_use ~ count +  GDD + precip + saturated_thickness | county_year + long_lat | 0 | long_lat, 
                             data = wells_within_buffer_radius_panel_4)

panel_FE_02_N_1_2_2016 <- felm(water_use ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                               data = wells_within_buffer_radius_panel_4)


panel_FE_02_avg_4 <- felm(water_use ~ year_2 +  GDD + precip + saturated_thickness + I(saturated_thickness^2) | county_year + long_lat | 0 | long_lat, 
                          data = wells_within_buffer_radius_panel_4)

panel_FE_02_N_1_4 <- felm(water_use ~ count +  GDD + precip + saturated_thickness + I(saturated_thickness^2) | county_year + long_lat | 0 | long_lat, 
                          data = wells_within_buffer_radius_panel_4)


###############################################################################################################################################
########################################################## Marginal Effects *9* mile ##########################################################
###############################################################################################################################################

foo = data.table()
qux = data.table()
wells_within_buffer_radius_panel_3 = data.table()

#-----------------
wells_within_buffer_radius_panel_3 = copy(wells_within_buffer_radius_panel_2[WUA_YEAR>2000])
setkey(wells_within_buffer_radius_panel_3, id)

qux = wells_within_buffer_radius_panel_3[,.(id, WUA_YEAR)]
foo = CREP_distance[distance <= 9 & distance > 0]
setkey(foo, id)
setkey(qux, id)
foo= foo[qux]
foo= foo[complete.cases(group)]
foo[, distance := round(distance, 3)]
foo[, count:= ifelse(signed_year_2+1 <= WUA_YEAR, 1, 0)]
foo[, min_distance  := min(distance), by=c("id", "WUA_YEAR")]
foo[, count := sum(count), by=c("id", "WUA_YEAR")]
foo = foo[distance == min_distance]
foo = unique(foo, by=c("id", "WUA_YEAR"))
foo = foo[,.(id, group, distance, count)]

setkey(foo, id)
setkey(wells_within_buffer_radius_panel_2, id)
wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_2[foo]

# wells_within_buffer_radius_panel_3 = unique(wells_within_buffer_radius_panel_3, by=c("PDIV_ID", "WUA_YEAR"))

wells_within_buffer_radius_panel_3[, group_year := paste(group, WUA_YEAR, sep = "-")]

wells_within_buffer_radius_panel_3[, long_lat := paste(LONGITUDE, LATITUDE, sep = "-")]
wells_within_buffer_radius_panel_3 = unique(wells_within_buffer_radius_panel_3, by=c("long_lat", "WUA_YEAR"))
wells_within_buffer_radius_panel_3[, year_2 := ifelse(count>0, 1, 0)]
wells_within_buffer_radius_panel_3[, dist_2 := 1/distance^4]

# saveRDS(wells_within_buffer_radius_panel_3, "regression_data/appendix_reg_table_1.rds")



panel_FE_09_count <- felm(water_use ~ count +  GDD + precip | long_lat + county_year | 0 | long_lat, weights = wells_within_buffer_radius_panel_3$dist_2, 
                          data = wells_within_buffer_radius_panel_3)

panel_FE_09_avg <- felm(water_use ~ year_2 +  GDD + precip | long_lat + county_year | 0 | long_lat, weights = wells_within_buffer_radius_panel_3$dist_2, 
                        data = wells_within_buffer_radius_panel_3)


###############################################################################################################################################
########################################################## Marginal Effects *3* mile ##########################################################
###############################################################################################################################################

foo = data.table()
qux = data.table()
wells_within_buffer_radius_panel_3 = data.table()

wells_within_buffer_radius_panel_3 = copy(wells_within_buffer_radius_panel_2[WUA_YEAR>2000])
setkey(wells_within_buffer_radius_panel_3, id)

qux = wells_within_buffer_radius_panel_3[,.(id, WUA_YEAR)]
foo = CREP_distance[distance <= 3 & distance > 0]
setkey(foo, id)
setkey(qux, id)
foo= foo[qux]
foo= foo[complete.cases(group)]
foo[, distance := round(distance, 3)]
foo[, count:= ifelse(signed_year_2+1 <= WUA_YEAR, 1, 0)]
foo[, min_distance  := min(distance), by=c("id", "WUA_YEAR")]
foo[, count := sum(count), by=c("id", "WUA_YEAR")]
foo = foo[distance == min_distance]
foo = unique(foo, by=c("id", "WUA_YEAR"))
foo = foo[,.(id, group, distance, count)]

setkey(foo, id)
setkey(wells_within_buffer_radius_panel_2, id)
wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_2[foo]


wells_within_buffer_radius_panel_3[, long_lat := paste(LONGITUDE, LATITUDE, sep = "-")]
wells_within_buffer_radius_panel_3 = unique(wells_within_buffer_radius_panel_3, by=c("long_lat", "WUA_YEAR"))

wells_within_buffer_radius_panel_3[, ACRES_IRR := as.numeric(ACRES_IRR)]
wells_within_buffer_radius_panel_3[, alfalfa := ifelse(CROP_CODE==1, 1, 0)]
wells_within_buffer_radius_panel_3[, corn    := ifelse(CROP_CODE==2, 1, 0)]
wells_within_buffer_radius_panel_3[, sorghum := ifelse(CROP_CODE==3, 1, 0)]
wells_within_buffer_radius_panel_3[, soybean := ifelse(CROP_CODE==4, 1, 0)]
wells_within_buffer_radius_panel_3[, wheat   := ifelse(CROP_CODE==5, 1, 0)]

wells_within_buffer_radius_panel_3[  !(range_K == "25 to 50" | range_K == "50 to 100" | range_K == "100 to 200") , range_K := NA]
wells_within_buffer_radius_panel_3[, conductivity := ifelse(is.na(range_K), NA, 
                                                            ifelse(range_K=="25 to 50" | range_K=="50 to 100", "low", "high"))]
wells_within_buffer_radius_panel_3[, conductivity :=  factor(conductivity)]
wells_within_buffer_radius_panel_3[, conductivity := relevel(conductivity, ref = "low")]


wells_within_buffer_radius_panel_3[, year_2 := ifelse(count>0, 1, 0)]
# 


wells_within_buffer_radius_panel_3[, year_2 := ifelse(count>0, 1, 0)]
# 

# saveRDS(wells_within_buffer_radius_panel_3, "regression_data/appendix_reg_table_1_2.rds")


panel_FE_03_avg <- felm(water_use ~ year_2 +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                        data = wells_within_buffer_radius_panel_3)

panel_FE_03_N_1 <- felm(water_use ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                        data = wells_within_buffer_radius_panel_3)


#=========================================================================================================================================================
#=========================                                          Diff in Diff                                          ================================
#=========================================================================================================================================================


###############################################################################################################################################
########################################################## Diff in Diff *2* mile ##########################################################
###############################################################################################################################################
foo = data.table()
qux = data.table()
wells_within_buffer_radius_panel_3 = data.table()
dt = data.table()

foo = copy(CREP_distance)
foo[, min_dist := min(distance), by="id"]
foo = foo[distance == min_dist]
foo[, InRing := ifelse(distance <= 2 & distance > 0, 1, 0)]
foo= unique(foo, by="id")
foo= foo[,.(id, InRing, CREP_well= group, signed_year_2, distance)]
foo = foo[distance < 15]

wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_2[,.(WR_ID, PDIV_ID, id, LONGITUDE, LATITUDE, WUA_YEAR, water_use, saturated_thickness, precip, GDD, county_year, TR_year, section_year_2, section_only, AUTH_QUANT, lower_K, SYSTEM, CROP_CODE, range_year)]

setkey(wells_within_buffer_radius_panel_3, id)
setkey(foo, id)
wells_within_buffer_radius_panel_3 = foo[wells_within_buffer_radius_panel_3]
wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_3[complete.cases(InRing)]
wells_within_buffer_radius_panel_3[, InRing_D   := InRing * distance]
wells_within_buffer_radius_panel_3[, PostRing   := ifelse(WUA_YEAR > signed_year_2, 1, 0)]
wells_within_buffer_radius_panel_3[, PostRing   := PostRing * InRing]
wells_within_buffer_radius_panel_3[, PostRing_D := PostRing * distance]
wells_within_buffer_radius_panel_3[, TPost      := WUA_YEAR - signed_year_2]
wells_within_buffer_radius_panel_3[, TPost      := ifelse(TPost < 0, 0, TPost)]
wells_within_buffer_radius_panel_3[, TPost      := TPost * InRing]
wells_within_buffer_radius_panel_3[, TPost_D    := TPost * distance]
wells_within_buffer_radius_panel_3[, sat_thickness_lag := c(NA,  saturated_thickness[-.N]), by=c("LONGITUDE", "LATITUDE")]
wells_within_buffer_radius_panel_3[, sat_thickness_diff:= saturated_thickness - sat_thickness_lag]
wells_within_buffer_radius_panel_3[, CREP_well  := CREP_well * InRing]
wells_within_buffer_radius_panel_3[, CREP_well  := as.character(CREP_well)]

# saveRDS(wells_within_buffer_radius_panel_3, "regression_data/appendix_reg_table_8.rds")

reg_diffindiff_2_mi = felm(water_use ~ PostRing + precip + GDD | factor(WUA_YEAR) + section_only + CREP_well | 0 | CREP_well,
                           data = wells_within_buffer_radius_panel_3)

reg_diffindiff_2_mi_b = felm(water_use ~ PostRing + precip + GDD  + AUTH_QUANT + lower_K + factor(SYSTEM) + CROP_CODE + LONGITUDE + LATITUDE | factor(WUA_YEAR) + section_only + CREP_well | 0 | CREP_well,
                             data = wells_within_buffer_radius_panel_3)

# ###############################################################################################################################################
# ########################################################## Diff in Diff *3* mile ##########################################################
# ###############################################################################################################################################
foo = data.table()
qux = data.table()
wells_within_buffer_radius_panel_3 = data.table()
dt = data.table()

foo = copy(CREP_distance)
foo[, min_dist := min(distance), by="id"]
foo = foo[distance == min_dist]
foo[, InRing := ifelse(distance <= 3 & distance > 0, 1, 0)]
foo= unique(foo, by="id")
foo= foo[,.(id, InRing, CREP_well= group, signed_year_2, distance)]
foo = foo[distance < 15]

wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_2[,.(WR_ID, PDIV_ID, id, LONGITUDE, LATITUDE, WUA_YEAR, water_use, saturated_thickness, precip, GDD, county_year, TR_year, section_year_2, section_only)]

setkey(wells_within_buffer_radius_panel_3, id)
setkey(foo, id)
wells_within_buffer_radius_panel_3 = foo[wells_within_buffer_radius_panel_3]
wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_3[complete.cases(InRing)]
wells_within_buffer_radius_panel_3[, InRing_D   := InRing * distance]
wells_within_buffer_radius_panel_3[, PostRing   := ifelse(WUA_YEAR > signed_year_2, 1, 0)]
wells_within_buffer_radius_panel_3[, PostRing   := PostRing * InRing]
wells_within_buffer_radius_panel_3[, PostRing_D := PostRing * distance]
wells_within_buffer_radius_panel_3[, TPost      := WUA_YEAR - signed_year_2]
wells_within_buffer_radius_panel_3[, TPost      := ifelse(TPost < 0, 0, TPost)]
wells_within_buffer_radius_panel_3[, TPost      := TPost * InRing]
wells_within_buffer_radius_panel_3[, TPost_D    := TPost * distance]
wells_within_buffer_radius_panel_3[, sat_thickness_lag := c(NA,  saturated_thickness[-.N]), by=c("LONGITUDE", "LATITUDE")]
wells_within_buffer_radius_panel_3[, sat_thickness_diff:= saturated_thickness - sat_thickness_lag]
wells_within_buffer_radius_panel_3[, CREP_well  := CREP_well * InRing]
wells_within_buffer_radius_panel_3[, CREP_well  := as.character(CREP_well)]

# saveRDS(wells_within_buffer_radius_panel_3, "regression_data/appendix_reg_table_8_2.rds")

reg_diffindiff_3_mi = felm(water_use ~ PostRing + precip + GDD | factor(WUA_YEAR) + section_only + CREP_well | 0 | CREP_well,
                           data = wells_within_buffer_radius_panel_3)


#=========================================================================================================================================================
#=========================                                             Non-pecuniary                                      ================================
#=========================================================================================================================================================
foo = data.table()
qux = data.table()
wells_within_buffer_radius_panel_3 = data.table()
dt = data.table()


#----------------- 3 miles
wells_within_buffer_radius_panel_3 = copy(wells_within_buffer_radius_panel_2[WUA_YEAR>2000])
setkey(wells_within_buffer_radius_panel_3, id)

qux = wells_within_buffer_radius_panel_3[,.(id, WUA_YEAR, WUAPERS_ID)]
foo = CREP_distance[distance <= 2 & distance > 0]
setkey(foo, id)
setkey(qux, id)
foo= foo[qux]
foo= foo[complete.cases(group)]
foo[, distance := round(distance, 3)]
foo[, count:= ifelse(signed_year_2+1 <= WUA_YEAR, 1, 0)]
foo[, min_distance  := min(distance), by=c("id", "WUA_YEAR")]
foo[, count := sum(count), by=c("id", "WUA_YEAR")]
foo[, AF_USED_2 := ifelse(signed_year_2+1 <= WUA_YEAR, AF_USED, 0)]
foo[, AF_USED_2 := sum(AF_USED_2), by=c("id", "WUA_YEAR")]
foo[, signed_year_2 := min(signed_year_2), by=c("WUAPERS_ID")]

foo[, count := mean(count), by=c("WUAPERS_ID", "WUA_YEAR")]
foo[, AF_USED_2 := mean(AF_USED_2), by=c("WUAPERS_ID", "WUA_YEAR")]

foo = unique(foo, by=c("WUAPERS_ID", "WUA_YEAR"))
foo = foo[,.(WUAPERS_ID, id_person = 1, WUA_YEAR, signed_year_2, count, AF_USED_2)]

setkey(foo, WUAPERS_ID, WUA_YEAR)
setkey(wells_within_buffer_radius_panel_3, WUAPERS_ID, WUA_YEAR)
wells_within_buffer_radius_panel_3 = foo[wells_within_buffer_radius_panel_3]

foo = copy(CREP_distance)
foo[, distance := min(distance), by="id"]
foo = foo[distance > 5 & distance < 15 ]
foo = unique(foo, by=c("id"))
foo = foo[,.(id, distance)]

setkey(foo, id)
setkey(wells_within_buffer_radius_panel_3, id)
wells_within_buffer_radius_panel_3 = foo[wells_within_buffer_radius_panel_3]
wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_3[complete.cases(county_year)]
wells_within_buffer_radius_panel_3[is.na(id_person), id_person := 0]

wells_within_buffer_radius_panel_3[, id_person := max(id_person), by="WUAPERS_ID"]
wells_within_buffer_radius_panel_3[,.N,by=id_person]

# wells_within_buffer_radius_panel_3[id_person==1 & is.na(count)]

wells_within_buffer_radius_panel_3[is.na(count), count := 0]
wells_within_buffer_radius_panel_3[, count := max(count), by=c("WUAPERS_ID", "WUA_YEAR")]
wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_3[id_person==1]

wells_within_buffer_radius_panel_3[, InRing := id_person]
wells_within_buffer_radius_panel_3[, PostRing   := ifelse(WUA_YEAR > 2008, 1, 0)]

wells_within_buffer_radius_panel_3[, year_2   := ifelse(count > 0, 1, 0)]
wells_within_buffer_radius_panel_3[, long_lat := paste(LONGITUDE, LATITUDE, sep = "-")]

# wells_within_buffer_radius_panel_3[, hist(distance)]

# saveRDS(wells_within_buffer_radius_panel_3, "regression_data/main_reg_table_3.rds")

panel_FE_02_avg_placebo_1 <- felm(water_use ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat,
                                  data = wells_within_buffer_radius_panel_3[distance > 4 & distance < 10])

panel_FE_02_avg_placebo_2 <- felm(water_use ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat,
                                  data = wells_within_buffer_radius_panel_3[distance > 10 & distance < 15])



#======================================================================
#==================       Event Study (with N)      ===================
#======================================================================

wells_within_buffer_radius_panel_3 = copy(wells_within_buffer_radius_panel_2[WUA_YEAR>2000])
setkey(wells_within_buffer_radius_panel_3, id)

qux = wells_within_buffer_radius_panel_3[,.(id, WUA_YEAR)]
foo = CREP_distance[distance <= 2 & distance > 0]
setkey(foo, id)
setkey(qux, id)
foo= foo[qux]
foo= foo[complete.cases(group)]
foo[, distance := round(distance, 3)]
foo[, total_wells:= ifelse(signed_year_2+1 <= WUA_YEAR, 1, 0)]
foo[, min_distance  := min(distance), by=c("id", "WUA_YEAR")]
foo[, total_wells := sum(total_wells), by=c("id", "WUA_YEAR")]
foo[, signed_year_2 := min(signed_year_2), by=c("id", "WUA_YEAR")]
foo[, year_2 := WUA_YEAR - signed_year_2]
foo = foo[distance == min_distance]
foo = unique(foo, by=c("id", "WUA_YEAR"))
foo = foo[,.(id, group, distance, signed_year_2, total_wells, year_2)]

setkey(foo, id)
setkey(wells_within_buffer_radius_panel_2, id)
wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_2[foo]

wells_within_buffer_radius_panel_3[, time_trend_event := year_2]
wells_within_buffer_radius_panel_3[time_trend_event < 0 , time_trend_event := 0]

# saveRDS(wells_within_buffer_radius_panel_3, "regression_data/main_reg_fig_5.rds")

wells_within_buffer_radius_panel_3[, time_trend_event := factor(time_trend_event)]
wells_within_buffer_radius_panel_3[, time_trend_event := relevel(time_trend_event, ref = "0")]
wells_within_buffer_radius_panel_3[, group_year := paste(group, WUA_YEAR, sep = "-")]
wells_within_buffer_radius_panel_3[, long_lat := paste(LONGITUDE, LATITUDE, sep = "-")]
wells_within_buffer_radius_panel_3 = unique(wells_within_buffer_radius_panel_3, by=c("long_lat", "WUA_YEAR"))

wells_within_buffer_radius_panel_3[, sat_thickness_lag := c(NA,  saturated_thickness[-.N]), by=c("LONGITUDE", "LATITUDE")]
wells_within_buffer_radius_panel_3[, sat_thickness_diff:= saturated_thickness - sat_thickness_lag]


mean_wells = round(wells_within_buffer_radius_panel_3[total_wells>0,mean(total_wells)], 2)

reg_1_ev = plm(water_use ~ total_wells + time_trend_event + time_trend_event:total_wells + GDD + precip + factor(county_year),
               index = c("long_lat", "WUA_YEAR"),
               data = wells_within_buffer_radius_panel_3, model = "within", effect="individual")

a_01_ev =summary(glht(reg_1_ev, linfct = c("total_wells * 2.85 + time_trend_event1 + total_wells:time_trend_event1  * 2.85 = 0")))[9]
a_02_ev =summary(glht(reg_1_ev, linfct = c("total_wells * 2.85 + time_trend_event2 + total_wells:time_trend_event2  * 2.85 = 0")))
a_03_ev =summary(glht(reg_1_ev, linfct = c("total_wells * 2.85 + time_trend_event3 + total_wells:time_trend_event3  * 2.85 = 0")))
a_04_ev =summary(glht(reg_1_ev, linfct = c("total_wells * 2.85 + time_trend_event4 + total_wells:time_trend_event4  * 2.85 = 0")))
a_05_ev =summary(glht(reg_1_ev, linfct = c("total_wells * 2.85 + time_trend_event5 + total_wells:time_trend_event5  * 2.85 = 0")))
a_06_ev =summary(glht(reg_1_ev, linfct = c("total_wells * 2.85 + time_trend_event6 + total_wells:time_trend_event6  * 2.85 = 0")))
a_07_ev =summary(glht(reg_1_ev, linfct = c("total_wells * 2.85 + time_trend_event7 + total_wells:time_trend_event7  * 2.85 = 0")))
a_08_ev =summary(glht(reg_1_ev, linfct = c("total_wells * 2.85 + time_trend_event8 + total_wells:time_trend_event8  * 2.85 = 0")))

fig1_ev = data.table(Estimate = c(a_01_ev$test$coefficients, 
                                  a_02_ev$test$coefficients, 
                                  a_03_ev$test$coefficients, 
                                  a_04_ev$test$coefficients, 
                                  a_05_ev$test$coefficients, 
                                  a_06_ev$test$coefficients, 
                                  a_07_ev$test$coefficients, 
                                  a_08_ev$test$coefficients
), 
std_err = c(a_01_ev$test$sigma, 
            a_02_ev$test$sigma, 
            a_03_ev$test$sigma, 
            a_04_ev$test$sigma, 
            a_05_ev$test$sigma, 
            a_06_ev$test$sigma, 
            a_07_ev$test$sigma, 
            a_08_ev$test$sigma
))

fig1_ev[, year := 1:8]
fig_point_estimates_ev=ggplot(fig1_ev, aes(x=factor(year), y=Estimate))+
  geom_point( size=1.5)+ 
  geom_errorbar(aes(ymin=Estimate-1.96*std_err, ymax=Estimate+1.96*std_err), width=.2,
                position=position_dodge(0.05))+
  geom_hline(yintercept = 0)+
  theme_bw()+
  xlab("Year")+
  ylab("Change in water use (AF)") +
  theme(axis.text    = element_text(size=12),
        axis.title   = element_text(size=12),
        legend.position   = "bottom",
        legend.text       = element_text(size=12),
        legend.title      = element_text(size=12),
        legend.key.size   = unit(1,"cm")
  )+ guides(col=guide_legend(title=""))

#=========================================================================================
#=========================       Panel FE for all wells   ================================
#=========================================================================================
# these are generated the same as above, except that they include ALL the wells


wells_within_buffer_radius_panel_2_all[, range_year   := paste(TWP, TWP_DIR, RNG, RNG_DIR, WUA_YEAR, sep = "-")]

foo = data.table()
qux = data.table()
wells_within_buffer_radius_panel_3 = data.table()

#-----------------
wells_within_buffer_radius_panel_3 = copy(wells_within_buffer_radius_panel_2_all[WUA_YEAR>2000])
setkey(wells_within_buffer_radius_panel_3, id)

qux = wells_within_buffer_radius_panel_3[,.(id, WUA_YEAR)]
foo = CREP_distance[distance <= 2 & distance > 0]
setkey(foo, id)
setkey(qux, id)
foo= foo[qux]
foo= foo[complete.cases(group)]
foo[, distance := round(distance, 3)]
foo[, count:= ifelse(signed_year_2+1 <= WUA_YEAR, 1, 0)]
foo[, min_distance  := min(distance), by=c("id", "WUA_YEAR")]
foo[, count := sum(count), by=c("id", "WUA_YEAR")]
foo = foo[distance == min_distance]
foo = unique(foo, by=c("id", "WUA_YEAR"))
foo = foo[,.(id, group, distance, count)]

setkey(foo, id)
setkey(wells_within_buffer_radius_panel_2_all, id)
wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_2_all[foo]
# 
wells_within_buffer_radius_panel_3[, group_year := paste(group, WUA_YEAR, sep = "-")]
wells_within_buffer_radius_panel_3[, long_lat := paste(LONGITUDE, LATITUDE, sep = "-")]
wells_within_buffer_radius_panel_3 = unique(wells_within_buffer_radius_panel_3, by=c("long_lat", "WUA_YEAR"))
wells_within_buffer_radius_panel_3[, ACRES_IRR := as.numeric(ACRES_IRR)]
wells_within_buffer_radius_panel_3[, alfalfa := ifelse(CROP_CODE==1, 1, 0)]
wells_within_buffer_radius_panel_3[, corn    := ifelse(CROP_CODE==2, 1, 0)]
wells_within_buffer_radius_panel_3[, sorghum := ifelse(CROP_CODE==3, 1, 0)]
wells_within_buffer_radius_panel_3[, soybean := ifelse(CROP_CODE==4, 1, 0)]
wells_within_buffer_radius_panel_3[, wheat   := ifelse(CROP_CODE==5, 1, 0)]
wells_within_buffer_radius_panel_3[  !(range_K == "25 to 50" | range_K == "50 to 100" | range_K == "100 to 200") , range_K := NA]
wells_within_buffer_radius_panel_3[, conductivity := ifelse(is.na(range_K), NA, 
                                                            ifelse(range_K=="25 to 50" | range_K=="50 to 100", "low", "high"))]
wells_within_buffer_radius_panel_3[, conductivity :=  factor(conductivity)]
wells_within_buffer_radius_panel_3[, conductivity := relevel(conductivity, ref = "low")]
wells_within_buffer_radius_panel_3[, year_2 := ifelse(count>0, 1, 0)]
# 

# saveRDS(wells_within_buffer_radius_panel_3, "regression_data/appendix_reg_table_1_3.rds")

panel_FE_02_avg_all <- felm(water_use ~ year_2 +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                            data = wells_within_buffer_radius_panel_3)

panel_FE_02_N_1_all <- felm(water_use ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                            data = wells_within_buffer_radius_panel_3)


#=========================================================================================
#===================================        Appendix      ================================
#=========================================================================================

# Placebo test: diff in diff
foo = data.table()
qux = data.table()
wells_within_buffer_radius_panel_3 = data.table()
dt = data.table()

foo = copy(CREP_distance)
foo[, min_dist := min(distance), by="id"]
foo = foo[distance == min_dist]
foo[, InRing := ifelse(distance <= 2 & distance > 0, 1, 0)]
foo= unique(foo, by="id")
foo= foo[,.(id, InRing, CREP_well= group, signed_year_2, distance)]
foo = foo[distance < 4]

wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_2[,.(WR_ID, PDIV_ID, id, LONGITUDE, LATITUDE, WUA_YEAR, water_use, saturated_thickness, precip, GDD, county_year, TR_year, section_year_2, section_only, AUTH_QUANT, lower_K, SYSTEM, CROP_CODE, range_year)]

setkey(wells_within_buffer_radius_panel_3, id)
setkey(foo, id)
wells_within_buffer_radius_panel_3 = foo[wells_within_buffer_radius_panel_3]
wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_3[complete.cases(InRing)]
wells_within_buffer_radius_panel_3[, InRing_D   := InRing * distance]
wells_within_buffer_radius_panel_3[, CREP_well  := CREP_well * InRing]
wells_within_buffer_radius_panel_3[, CREP_well  := as.character(CREP_well)]
wells_within_buffer_radius_panel_3[, treat_year := ifelse(WUA_YEAR > 2001, 1, 0)]

# saveRDS(wells_within_buffer_radius_panel_3, "regression_data/appendix_fig_1.rds")

reg_diffindiff_2_mi_2 = felm(water_use ~ InRing:factor(WUA_YEAR) + precip + GDD + saturated_thickness | factor(WUA_YEAR) + section_only + CREP_well | 0 | CREP_well,
                             data = wells_within_buffer_radius_panel_3[WUA_YEAR < 2008])

coefs_DiD = data.table(year=2002:2007, coefficient_DiD = rep(0 , 6), standard_error_DiD= rep(0, 6))
b = data.table()

for(i in 2002:2007){
  reg_diffindiff_2_mi_2 = felm(water_use ~ InRing + treat_year + InRing:treat_year + precip + GDD + saturated_thickness,
                               data = wells_within_buffer_radius_panel_3[WUA_YEAR == 2001 | WUA_YEAR == i])
  beta_DiD = summary(reg_diffindiff_2_mi_2)$coefficients[nrow(summary(reg_diffindiff_2_mi_2)$coefficients),1]
  std_error= summary(reg_diffindiff_2_mi_2)$coefficients[nrow(summary(reg_diffindiff_2_mi_2)$coefficients),2]
  coefs_DiD[i-2001, 2] = beta_DiD
  coefs_DiD[i-2001, 3] = std_error
  # print(i)
}

fig_placebo_diff_in_diff = ggplot(data = coefs_DiD, aes(x=factor(year), y=coefficient_DiD, group=1))+
  geom_errorbar(aes(ymin=coefficient_DiD - 1.96 * standard_error_DiD , ymax=coefficient_DiD + 1.96 * standard_error_DiD), width=.1) +
  geom_line() +
  geom_point() +
  theme_bw()+
  # geom_errorbar(data = dt_ESE, aes(y=mean_extra_discharge, ymin=mean_extra_discharge - 1.96 * sd_extra_discharge , ymax=mean_extra_discharge + 1.96 * sd_extra_discharge),
  #               width=.2, col="red") +
  xlab("Year")+
  ylab("Diff in Diff Coefficient")+
  theme(axis.text    = element_text(size=24),
        axis.title   = element_text(size=24),
        axis.text.x  = element_text(angle = 90, hjust = 1),
        legend.position   = "bottom",
        legend.text       = element_text(size=24),
        legend.title      = element_text(size=24),
        legend.key.size   = unit(1,"cm")
  )  






# Spatial dissipation of effects

# regressions with different mile bands
# ---
# -   -   -  2 to 5
# ---
foo = data.table()
qux = data.table()
wells_within_buffer_radius_panel_3 = data.table()

wells_within_buffer_radius_panel_3 = copy(wells_within_buffer_radius_panel_2[WUA_YEAR>2000])
setkey(wells_within_buffer_radius_panel_3, id)

qux = wells_within_buffer_radius_panel_3[,.(id, WUA_YEAR)]
foo = CREP_distance[distance <= 5 & distance > 2]
setkey(foo, id)
setkey(qux, id)
foo= foo[qux]
foo= foo[complete.cases(group)]
foo[, distance := round(distance, 3)]
foo[, count:= ifelse(signed_year_2+1 <= WUA_YEAR, 1, 0)]
foo[, min_distance  := min(distance), by=c("id", "WUA_YEAR")]
foo[, count := sum(count), by=c("id", "WUA_YEAR")]
foo = foo[distance == min_distance]
foo = unique(foo, by=c("id", "WUA_YEAR"))
foo = foo[,.(id, group, distance, count)]

setkey(foo, id)
setkey(wells_within_buffer_radius_panel_2, id)
wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_2[foo]

# wells_within_buffer_radius_panel_3 = unique(wells_within_buffer_radius_panel_3, by=c("PDIV_ID", "WUA_YEAR"))

wells_within_buffer_radius_panel_3[, long_lat := paste(LONGITUDE, LATITUDE, sep = "-")]
wells_within_buffer_radius_panel_3 = unique(wells_within_buffer_radius_panel_3, by=c("long_lat", "WUA_YEAR"))

wells_within_buffer_radius_panel_3[, ACRES_IRR := as.numeric(ACRES_IRR)]
wells_within_buffer_radius_panel_3[, alfalfa := ifelse(CROP_CODE==1, 1, 0)]
wells_within_buffer_radius_panel_3[, corn    := ifelse(CROP_CODE==2, 1, 0)]
wells_within_buffer_radius_panel_3[, sorghum := ifelse(CROP_CODE==3, 1, 0)]
wells_within_buffer_radius_panel_3[, soybean := ifelse(CROP_CODE==4, 1, 0)]
wells_within_buffer_radius_panel_3[, wheat   := ifelse(CROP_CODE==5, 1, 0)]

wells_within_buffer_radius_panel_3[  !(range_K == "25 to 50" | range_K == "50 to 100" | range_K == "100 to 200") , range_K := NA]
wells_within_buffer_radius_panel_3[, conductivity := ifelse(is.na(range_K), NA, 
                                                            ifelse(range_K=="25 to 50" | range_K=="50 to 100", "low", "high"))]
wells_within_buffer_radius_panel_3[, conductivity :=  factor(conductivity)]
wells_within_buffer_radius_panel_3[, conductivity := relevel(conductivity, ref = "low")]
wells_within_buffer_radius_panel_3[, year_2 := ifelse(count>0, 1, 0)]

# saveRDS(wells_within_buffer_radius_panel_3, "regression_data/appendix_table_2_1.rds")

panel_FE_03_avg_2_5 <- felm(water_use ~ year_2 +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                            data = wells_within_buffer_radius_panel_3)

panel_FE_03_N_1_2_5 <- felm(water_use ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                            data = wells_within_buffer_radius_panel_3)

# ---
# -   -   -  5 to 8
# ---
foo = data.table()
qux = data.table()
wells_within_buffer_radius_panel_3 = data.table()

wells_within_buffer_radius_panel_3 = copy(wells_within_buffer_radius_panel_2[WUA_YEAR>2000])
setkey(wells_within_buffer_radius_panel_3, id)

qux = wells_within_buffer_radius_panel_3[,.(id, WUA_YEAR)]
foo = CREP_distance[distance <= 8 & distance > 5]
setkey(foo, id)
setkey(qux, id)
foo= foo[qux]
foo= foo[complete.cases(group)]
foo[, distance := round(distance, 3)]
foo[, count:= ifelse(signed_year_2+1 <= WUA_YEAR, 1, 0)]
foo[, min_distance  := min(distance), by=c("id", "WUA_YEAR")]
foo[, count := sum(count), by=c("id", "WUA_YEAR")]
foo = foo[distance == min_distance]
foo = unique(foo, by=c("id", "WUA_YEAR"))
foo = foo[,.(id, group, distance, count)]

setkey(foo, id)
setkey(wells_within_buffer_radius_panel_2, id)
wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_2[foo]

# wells_within_buffer_radius_panel_3 = unique(wells_within_buffer_radius_panel_3, by=c("PDIV_ID", "WUA_YEAR"))

wells_within_buffer_radius_panel_3[, long_lat := paste(LONGITUDE, LATITUDE, sep = "-")]
wells_within_buffer_radius_panel_3 = unique(wells_within_buffer_radius_panel_3, by=c("long_lat", "WUA_YEAR"))

wells_within_buffer_radius_panel_3[, ACRES_IRR := as.numeric(ACRES_IRR)]
wells_within_buffer_radius_panel_3[, alfalfa := ifelse(CROP_CODE==1, 1, 0)]
wells_within_buffer_radius_panel_3[, corn    := ifelse(CROP_CODE==2, 1, 0)]
wells_within_buffer_radius_panel_3[, sorghum := ifelse(CROP_CODE==3, 1, 0)]
wells_within_buffer_radius_panel_3[, soybean := ifelse(CROP_CODE==4, 1, 0)]
wells_within_buffer_radius_panel_3[, wheat   := ifelse(CROP_CODE==5, 1, 0)]

wells_within_buffer_radius_panel_3[  !(range_K == "25 to 50" | range_K == "50 to 100" | range_K == "100 to 200") , range_K := NA]
wells_within_buffer_radius_panel_3[, conductivity := ifelse(is.na(range_K), NA, 
                                                            ifelse(range_K=="25 to 50" | range_K=="50 to 100", "low", "high"))]
wells_within_buffer_radius_panel_3[, conductivity :=  factor(conductivity)]
wells_within_buffer_radius_panel_3[, conductivity := relevel(conductivity, ref = "low")]
wells_within_buffer_radius_panel_3[, year_2 := ifelse(count>0, 1, 0)]


# saveRDS(wells_within_buffer_radius_panel_3, "regression_data/appendix_table_2_2.rds")


panel_FE_03_avg_5_8 <- felm(water_use ~ year_2 +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                            data = wells_within_buffer_radius_panel_3)

panel_FE_03_N_1_5_8 <- felm(water_use ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                            data = wells_within_buffer_radius_panel_3)


# ---
# -   -   -  8 to 11
# ---
foo = data.table()
qux = data.table()
wells_within_buffer_radius_panel_3 = data.table()

wells_within_buffer_radius_panel_3 = copy(wells_within_buffer_radius_panel_2[WUA_YEAR>2000])
setkey(wells_within_buffer_radius_panel_3, id)

qux = wells_within_buffer_radius_panel_3[,.(id, WUA_YEAR)]
foo = CREP_distance[distance <= 11 & distance > 8]
setkey(foo, id)
setkey(qux, id)
foo= foo[qux]
foo= foo[complete.cases(group)]
foo[, distance := round(distance, 3)]
foo[, count:= ifelse(signed_year_2+1 <= WUA_YEAR, 1, 0)]
foo[, min_distance  := min(distance), by=c("id", "WUA_YEAR")]
foo[, count := sum(count), by=c("id", "WUA_YEAR")]
foo = foo[distance == min_distance]
foo = unique(foo, by=c("id", "WUA_YEAR"))
foo = foo[,.(id, group, distance, count)]

setkey(foo, id)
setkey(wells_within_buffer_radius_panel_2, id)
wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_2[foo]

# wells_within_buffer_radius_panel_3 = unique(wells_within_buffer_radius_panel_3, by=c("PDIV_ID", "WUA_YEAR"))

wells_within_buffer_radius_panel_3[, long_lat := paste(LONGITUDE, LATITUDE, sep = "-")]
wells_within_buffer_radius_panel_3 = unique(wells_within_buffer_radius_panel_3, by=c("long_lat", "WUA_YEAR"))

wells_within_buffer_radius_panel_3[, ACRES_IRR := as.numeric(ACRES_IRR)]
wells_within_buffer_radius_panel_3[, alfalfa := ifelse(CROP_CODE==1, 1, 0)]
wells_within_buffer_radius_panel_3[, corn    := ifelse(CROP_CODE==2, 1, 0)]
wells_within_buffer_radius_panel_3[, sorghum := ifelse(CROP_CODE==3, 1, 0)]
wells_within_buffer_radius_panel_3[, soybean := ifelse(CROP_CODE==4, 1, 0)]
wells_within_buffer_radius_panel_3[, wheat   := ifelse(CROP_CODE==5, 1, 0)]

wells_within_buffer_radius_panel_3[  !(range_K == "25 to 50" | range_K == "50 to 100" | range_K == "100 to 200") , range_K := NA]
wells_within_buffer_radius_panel_3[, conductivity := ifelse(is.na(range_K), NA, 
                                                            ifelse(range_K=="25 to 50" | range_K=="50 to 100", "low", "high"))]
wells_within_buffer_radius_panel_3[, conductivity :=  factor(conductivity)]
wells_within_buffer_radius_panel_3[, conductivity := relevel(conductivity, ref = "low")]
wells_within_buffer_radius_panel_3[, year_2 := ifelse(count>0, 1, 0)]


# saveRDS(wells_within_buffer_radius_panel_3, "regression_data/appendix_table_2_3.rds")


panel_FE_03_avg_8_11 <- felm(water_use ~ year_2 +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                             data = wells_within_buffer_radius_panel_3)

panel_FE_03_N_1_8_11 <- felm(water_use ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                             data = wells_within_buffer_radius_panel_3)

# ---
# -   -   -  11 - 14
# ---
foo = data.table()
qux = data.table()
wells_within_buffer_radius_panel_3 = data.table()

wells_within_buffer_radius_panel_3 = copy(wells_within_buffer_radius_panel_2[WUA_YEAR>2000])
setkey(wells_within_buffer_radius_panel_3, id)

qux = wells_within_buffer_radius_panel_3[,.(id, WUA_YEAR)]
foo = CREP_distance[distance <= 14 & distance > 11]
setkey(foo, id)
setkey(qux, id)
foo= foo[qux]
foo= foo[complete.cases(group)]
foo[, distance := round(distance, 3)]
foo[, count:= ifelse(signed_year_2+1 <= WUA_YEAR, 1, 0)]
foo[, min_distance  := min(distance), by=c("id", "WUA_YEAR")]
foo[, count := sum(count), by=c("id", "WUA_YEAR")]
foo = foo[distance == min_distance]
foo = unique(foo, by=c("id", "WUA_YEAR"))
foo = foo[,.(id, group, distance, count)]

setkey(foo, id)
setkey(wells_within_buffer_radius_panel_2, id)
wells_within_buffer_radius_panel_3 = wells_within_buffer_radius_panel_2[foo]

# wells_within_buffer_radius_panel_3 = unique(wells_within_buffer_radius_panel_3, by=c("PDIV_ID", "WUA_YEAR"))

wells_within_buffer_radius_panel_3[, long_lat := paste(LONGITUDE, LATITUDE, sep = "-")]
wells_within_buffer_radius_panel_3 = unique(wells_within_buffer_radius_panel_3, by=c("long_lat", "WUA_YEAR"))

wells_within_buffer_radius_panel_3[, ACRES_IRR := as.numeric(ACRES_IRR)]
wells_within_buffer_radius_panel_3[, alfalfa := ifelse(CROP_CODE==1, 1, 0)]
wells_within_buffer_radius_panel_3[, corn    := ifelse(CROP_CODE==2, 1, 0)]
wells_within_buffer_radius_panel_3[, sorghum := ifelse(CROP_CODE==3, 1, 0)]
wells_within_buffer_radius_panel_3[, soybean := ifelse(CROP_CODE==4, 1, 0)]
wells_within_buffer_radius_panel_3[, wheat   := ifelse(CROP_CODE==5, 1, 0)]

wells_within_buffer_radius_panel_3[  !(range_K == "25 to 50" | range_K == "50 to 100" | range_K == "100 to 200") , range_K := NA]
wells_within_buffer_radius_panel_3[, conductivity := ifelse(is.na(range_K), NA, 
                                                            ifelse(range_K=="25 to 50" | range_K=="50 to 100", "low", "high"))]
wells_within_buffer_radius_panel_3[, conductivity :=  factor(conductivity)]
wells_within_buffer_radius_panel_3[, conductivity := relevel(conductivity, ref = "low")]
wells_within_buffer_radius_panel_3[, year_2 := ifelse(count>0, 1, 0)]


# saveRDS(wells_within_buffer_radius_panel_3, "regression_data/appendix_table_2_4.rds")


panel_FE_03_avg_11_14 <- felm(water_use ~ year_2 +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                              data = wells_within_buffer_radius_panel_3)

panel_FE_03_N_1_11_14 <- felm(water_use ~ count +  GDD + precip | county_year + long_lat | 0 | long_lat, 
                              data = wells_within_buffer_radius_panel_3)



```






\newpage

# Figures and Tables

Please note that table 1 of the manuscript is not added here so that table numbers are off by 1, e.g., Table 1 here is Table 2 in the manuscript.


```{r, results='asis', echo=F}
stargazer(sum_stats, summary=F, rownames = F, header = F, title = "\\label{tab:summary_stats}Summary statistics for CREP wells, wells within a 2-mile radius of CREP wells and wells outside the 2-mile radius of CREP well before 2008.")


```

\newpage

\begin{landscape}


```{r, results='asis', echo=F}
stargazer(panel_FE_02_avg_2, panel_FE_02_avg_2_2016, panel_FE_02_avg_2016, panel_FE_02_avg_3, panel_FE_02_N_1_2,  panel_FE_02_N_1_2_2016, panel_FE_02_N_1_2016, panel_FE_02_N_1_3,  no.space=TRUE, header = F, 
          dep.var.caption = "Groundwater Extraction (AF)", dep.var.labels.include = FALSE,
          column.sep.width = ".1pt", 
          covariate.labels = c("Post treatment", "Number of CREP wells", "Growing degree days", "Growing season precipitation", "Saturated thickness"
          ),
          omit = c("county_year", "section_only"),
          add.lines=list(c("County by year fixed effects?", "Yes", "Yes", "Yes", "No", "Yes", "Yes", "Yes", "No")
                         , c("CREP well by year fixed effects?", "No", "No", "No", "Yes", "No", "No", "No", "Yes")
                         , c("Well fixed effects?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")
          ),
          title = "\\label{tab:main_result}The average (columns 1, 2, 3 and 4) and marginal (columns 5, 6, 7 and 8) effects of the CREP program on water use of neighboring active wells within a 2-mile radius of a CREP well.", 
          notes = "Standard errors are clustered at the individual level",
          omit.stat = c("f", "ser")
)


```


\newpage


```{r, results='asis', echo=F}
stargazer(panel_FE_02_N_1_ST, panel_FE_02_cond, panel_FE_02_cond_3, panel_FE_02_avg_placebo_1, panel_FE_02_avg_placebo_2, no.space=TRUE, header = F, 
          dep.var.caption = "", dep.var.labels.include = FALSE,
          column.sep.width = ".1pt", 
          covariate.labels = c("Number of CREP wells", "Post treatment", "Growing degree days", "Growing season precipitation", #"Saturated thickness", 
                               "Number of CREP wells $\\times$ high conductivity", "Post treatment $\\times$ high conductivity"
          ),
          omit = c("county_year", "section_only"),
          add.lines=list(c("County by year fixed effects?", "Yes", "Yes", "Yes", "Yes", "Yes")
                         , c("Well fixed effects?", "Yes", "Yes", "Yes", "Yes", "Yes")
          ),
          title = "\\label{tab:mech_CPR}The effect of the number of wells in CREP on saturated thickness of neighboring active wells (column 1) and the effect of hydraulic conductivity on groundwater use of neighboring active wells within a 2-mile radius of a CREP well (columns 2 and 3). Columns 4 and 5 show the changes in water use at wells owned by individuals who have at least one well within a 2-mile radius of a CREP well.",
          omit.stat = c("f", "ser"), 
          # model.numbers = F,
          column.labels = c("Saturated thickness", "Extraction", "Extraction", "Ext. 4-10mi", "Ext. 10-15mi")
)


```

\end{landscape}

\newpage

\begin{landscape}

```{r, results='asis', echo=F}
stargazer(panel_FE_02_N_1_ACRES_IRR, panel_FE_02_N_1_ACRES_IRR_2, panel_FE_02_N_1_alfalfa, panel_FE_02_N_1_corn,  panel_FE_02_N_1_sorghum,  panel_FE_02_N_1_soybean,  
          panel_FE_02_N_1_wheat,  
          no.space=TRUE, header = F, 
          dep.var.caption = c(NULL, "Probability of growing"), dep.var.labels.include = FALSE,
          column.labels = c("Irr acres", "GW application", "Alfalfa", "Corn", "Sorghum", "Soybean", "Wheat"),
          column.sep.width = ".1pt", 
          covariate.labels = c("Number of CREP wells", "Growing degree days", "Growing season precipitation", #"Saturated thickness", 
                               "Irrigated acres"
          ),
          omit = c("county_year", "section_only"),
          add.lines=list(c("County by year fixed effects?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")
                         , c("Well fixed effects?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")
          ),
          title = "\\label{tab:mech_margins}The effect of the number of wells in CREP on irrigated acres, intensive margin application and probability of growing each crop for neighboring active wells within a 2-mile radius.", 
          notes = "Standard errors are clustered at the individual level",
          omit.stat = c("f", "ser")
)


```

\end{landscape}

\newpage


```{r, figCREParea, fig.show="hold", echo=F}
fig_CREP_area
```

\newpage

```{r, figCREPwells, fig.show="hold", echo=F}
fig_CREP_wells
```

\newpage

```{r, echo=FALSE, warning=FALSE, fig.cap= "\\label{fig:hist}Cumulative number of wells and percent of total wells that have enrolled in CREP and ceased groundwater extraction since 2008.", fig.show='hold', fig.align='center'}
fig_hist

```

\newpage

```{r, echo=FALSE, warning=FALSE, fig.cap= "\\label{fig:diff_in_diff}Figure shows average annual groundwater use for wells that are within 1-, 2-, and 3-mile radii of a CREP well.",fig.show='hold',fig.align='center'}
fig_motivation

```


\newpage


```{r, echo=FALSE, warning=FALSE, fig.cap= "\\label{fig:N_years}Figure shows the joint estimates of the effects of number of wells enrolled in CREP and years post treatment on groundwater use." ,fig.show='hold',fig.align='center', out.width= "60%"}
fig_point_estimates_ev

```

\newpage



# (APPENDIX) Appendix {-}

# Appendix A

\begin{landscape}

```{r, results='asis', echo=F}
stargazer(panel_FE_03_avg, panel_FE_02_avg_all, panel_FE_09_avg, panel_FE_03_N_1, panel_FE_02_N_1_all,  panel_FE_09_count,  no.space=TRUE, header = F, 
          dep.var.caption = "Groundwater Extraction (AF)", dep.var.labels.include = FALSE,
          column.sep.width = ".1pt", 
          covariate.labels = c("Post treatment", "Number of CREP wells", "Growing degree days", "Growing season precipitation", "Saturated thickness"
          ),
          omit = c("county_year", "section_only"),
          add.lines=list(c("County by year fixed effects?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")
                         , c("Well fixed effects?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")
          ),
          title = "\\label{tab:robustness_spec}The average (columns 1, 2 and 3) and marginal (columns 4, 5 and 6) effects of the CREP program on water use of neighboring active wells within a 2-mile radius. Columns 1 and 4 show the effects for wells within a 3-mi radius. Columns 2 and 5 show the effects for wells within a 2-mile radius including all the wells that are owned by individuals who enroll a well in CREP. Columns 3 and 6 show the effects for all the wells within a 9-mile radius where all the observation are weighted by the inverse square of the distance.", 
          notes = "Standard errors are clustered at the individual level",
          omit.stat = c("f", "ser"),
          column.labels = c("3 mi", "CREP", "9 mi", "3 mi", "CREP", "9 mi")
          
)


```


\newpage


```{r, results='asis', echo=F}
stargazer(panel_FE_03_avg_2_5, panel_FE_03_avg_5_8, panel_FE_03_avg_8_11, panel_FE_03_avg_11_14,  
          panel_FE_03_N_1_2_5, panel_FE_03_N_1_5_8, panel_FE_03_N_1_8_11, panel_FE_03_N_1_11_14,
          no.space=TRUE, header = F, 
          dep.var.caption = "Groundwater Extraction (AF)", dep.var.labels.include = FALSE,
          column.sep.width = ".1pt", 
          column.labels = c("(2-5)", "(5-8)", "(8-11)", "(11-14)", "(2-5)", "(5-8)", "(8-11)", "(11-14)"), 
          covariate.labels = c("Post treatment", "Number of CREP wells", "Growing degree days", "Growing season precipitation"
          ),
          omit = c("county_year", "section_only"),
          add.lines=list(c("County by year fixed effects?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")
                         , c("Well fixed effects?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")
          ),
          title = "\\label{tab:robustness_distance}The average (columns 1 to 4) and marginal (columns 5 to 8) effects of the CREP program on water use of the neighboring active wells within 2-mile radius for different distance bands.", 
          # notes.align = "l", 
          notes = "Standard errors are clustered at the individual level",
          omit.stat = c("f", "ser")
)
```

\end{landscape}

\newpage


```{r, results='asis', echo=F}
stargazer(panel_FE_02_N_1_2_factor,  no.space=TRUE, header = F, 
          dep.var.caption = "Groundwater Extraction (AF)", dep.var.labels.include = FALSE,
          column.sep.width = ".1pt", 
          covariate.labels = c("1 CREP well", 
                               "2 CREP wells", 
                               "3 CREP wells", 
                               "4 CREP wells", 
                               "5 CREP wells", 
                               ">5 CREP wells", 
                               "Growing degree days", "Growing season precipitation"
          ),
          omit = c("county_year", "section_only"),
          add.lines=list(c("County by year fixed effects?", "Yes")
                         , c("Well fixed effects?", "Yes")
          ),
          title = "\\label{tab:main_result_count}Marginal effects of the CREP program on water use of the neighboring active wells within 2-mile radius.", 
          # notes.align = "l", 
          notes = "Standard errors are clustered at the individual level",
          omit.stat = c("f", "ser")
)
```



\newpage


```{r, results='asis', echo=F}
stargazer(panel_FE_02_avg_2_org, panel_FE_02_avg_3_org, panel_FE_02_N_1_2_org,  panel_FE_02_N_1_3_org,  no.space=TRUE, header = F, 
          dep.var.caption = "Groundwater Extraction (AF)", dep.var.labels.include = FALSE,
          column.sep.width = ".1pt", 
          covariate.labels = c("Post treatment", "Number of CREP wells", "Growing degree days", "Growing season precipitation"
          ),
          omit = c("county_year", "section_only"),
          add.lines=list(c("County by year fixed effects?", "Yes", "No", "Yes", "No")
                         , c("CREP well by year fixed effects?", "No", "Yes", "No", "Yes")
                         , c("Well fixed effects?", "Yes", "Yes", "Yes", "Yes")
          ),
          title = "\\label{tab:robustness_CREP_all}The average (columns 1 and 2) and marginal (columns 3 and 4) effects of the CREP program on water use of the neighboring active wells within 2-mile radius. The sample for this regression includes all the wells that are owned by producers who enroll at least one well in the CREP program.", 
          # notes.align = "l", 
          notes = "Standard errors are clustered at the individual level",
          omit.stat = c("f", "ser")
)
```


\newpage



```{r, results='asis', echo=F}
stargazer(panel_FE_02_avg_2_1, panel_FE_02_avg_3_1, panel_FE_02_N_1_2_1,  panel_FE_02_N_1_3_1,  no.space=TRUE, header = F, 
          dep.var.caption = "Groundwater Extraction (AF)", dep.var.labels.include = FALSE,
          column.sep.width = ".1pt", 
          covariate.labels = c("Post treatment", "Number of CREP wells", "Growing degree days", "Growing season precipitation"
          ),
          omit = c("county_year", "section_only"),
          add.lines=list(c("County by year fixed effects?", "Yes", "No", "Yes", "No")
                         , c("CREP well by year fixed effects?", "No", "Yes", "No", "Yes")
                         , c("Well fixed effects?", "Yes", "Yes", "Yes", "Yes")
          ),
          title = "\\label{tab:robustness_CREP_CREP}The average (columns 1 and 2) and marginal (columns 3 and 4) effects of the CREP program on water use of the neighboring active wells within 2-mile radius. The sample for this regression includes all the wells that are owned by producers who enroll at least one well within a 2-mile radius of wells in the CREP program.", 
          # notes.align = "l", 
          notes = "Standard errors are clustered at the individual level",
          omit.stat = c("f", "ser")
)

```



\newpage


```{r, results='asis', echo=F}
stargazer(panel_FE_02_avg_2_0, panel_FE_02_avg_3_0, panel_FE_02_N_1_2_0,  panel_FE_02_N_1_3_0,  no.space=TRUE, header = F, 
          dep.var.caption = "Groundwater Extraction (AF)", dep.var.labels.include = FALSE,
          column.sep.width = ".1pt", 
          covariate.labels = c("Post treatment", "Number of CREP wells", "Growing degree days", "Growing season precipitation"
          ),
          omit = c("county_year", "section_only"),
          add.lines=list(c("County by year fixed effects?", "Yes", "No", "Yes", "No")
                         , c("CREP well by year fixed effects?", "No", "Yes", "No", "Yes")
                         , c("Well fixed effects?", "Yes", "Yes", "Yes", "Yes")
          ),
          title = "\\label{tab:robustness_CREP_no_CREP}The average (columns 1 and 2) and marginal (columns 3 and 4) effects of the CREP program on water use of the neighboring active wells within 2-mile radius. The sample for this regression includes all the wells that are owned by producers who enroll at least one well outside a 2-mile radius of wells in the CREP program.", 
          # notes.align = "l", 
          notes = "Standard errors are clustered at the individual level",
          omit.stat = c("f", "ser")
)
```



\newpage


\begin{landscape}


```{r, results='asis', echo=F}
stargazer(panel_FE_02_avg_4, panel_FE_02_avg_5, panel_FE_02_avg_6, panel_FE_02_N_1_4,  panel_FE_02_N_1_5,  panel_FE_02_N_1_6,  no.space=TRUE, header = F, 
          dep.var.caption = "Groundwater Extraction (AF)", dep.var.labels.include = FALSE,
          column.sep.width = ".1pt", 
          covariate.labels = c("Post treatment", "AF retired", "Number of CREP wells", "Ratio of CREP wells", "Growing degree days", "Growing season precipitation", "Saturated thickness", "Saturated thickness squared"
          ),
          omit = c("county_year", "section_only"),
          add.lines=list(c("County by year fixed effects?", "Yes", "Yes", "No", "Yes", "Yes", "No")
                         , c("Well fixed effects?", "Yes", "No", "Yes", "Yes", "No", "Yes")
                         , c("Well id fixed effects?", "No", "Yes", "No", "No", "Yes", "No")
          ),
          title = "\\label{tab:robustness_spec_2}The average (columns 1 and 2) and marginal (columns 4 and 5) effects of the CREP program on water use of neighboring active wells within a 2-mile radius. Columns 1 and 4 show the effects of treatment including the square of the saturated thickness. Columns 2 and 5 include well identifier as the unit of observation rather than longitude and latitude. Column 3 shows the effect of an additional acre-foot retired. Column 6 shows the effect of the ratio of wells retired.", 
          notes = "Standard errors are clustered at the individual level",
          omit.stat = c("f", "ser")
)


```

\end{landscape}

\newpage


```{r, results='asis', echo=F}
stargazer(reg_diffindiff_2_mi, reg_diffindiff_2_mi_b, reg_diffindiff_3_mi,  no.space=TRUE, header = F, 
          dep.var.caption = "Groundwater Extraction (AF)", dep.var.labels.include = FALSE,
          column.sep.width = ".1pt", 
          covariate.labels = c("Inside buffer $\\times$ post treatment", "Growing season precipitation", "Growing degree days", "Saturated thickness"
          ),
          omit = c("WUA_YEAR", "section_only", "CREP_well", "AUTH_QUANT", "lower_K", "SYSTEM", "CROP_CODE", "LONGITUDE", "LATITUDE"),
          add.lines=list(c("Year fixed effects?",      "Yes", "Yes",   "Yes")
                         , c("Section fixed effects",    "Yes", "Yes",   "Yes")
                         , c("CREP well fixed effects",  "Yes", "Yes",   "Yes")
                         , c("Other controls", "", "Crop type"          ,  "")
                         , c(""              , "", "Technology"         ,  "")
                         , c(""              , "", "Conductivity"        ,  "")
                         , c(""              , "", "Authorized quantity",  "")
                         , c(""              , "", "Longitude, Latitude",  "")
          ),
          title = "\\label{tab:robustness_DID}Difference-in-differences estimates of the effect of CREP enrollment on water use of active wells that are within a 2-mile radius of CREP wells.", 
          notes = "Standard errors are clustered at the CREP well level",
          omit.stat = c("f", "ser"),
          column.labels = c("2 miles", "2 miles", "3 miles", "3 miles")
)


```

\newpage




```{r, echo=FALSE, warning=FALSE, fig.cap= " The placebo test showing the difference-in-differences coefficient for groundwater use of wells within 2 miles of CREP wells (treated) and those between 2 and 4 miles of CREP wells (control) before the treatment period, i.e., 2001-2008. Th pseudo treatment year is set at 2002. Control variables are precipitation, growing degree days and saturated thickness." ,fig.show='hold',fig.align='center'}
fig_placebo_diff_in_diff

```











